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ABSTRACT 

We present a series of hydrodynamic simulations of isolated galaxies with stellar 
mass of 10 9 M©. The models use a resolution of 750 M© per particle and include a 
treatment for the full non-equilibrium chemical evolution of ions and molecules (157 
species in total), along with gas cooling rates computed self-consistently using the 
non-equilibrium abundances. We compare these to simulations evolved using cooling 
rates calculated assuming chemical (including ionisation) equilibrium, and we consider 
a wide range of metallicities and UV radiation fields, including a local prescription for 
self-shielding by gas and dust. We find higher star formation rates and stronger out¬ 
flows at higher metallicity and for weaker radiation fields, as gas can more easily cool 
to a cold (few hundred Kelvin) star forming phase under such conditions. Contrary 
to variations in the metallicity and the radiation field, non-equilibrium chemistry gen¬ 
erally has no strong effect on the total star formation rates or outflow properties. 
However, it is important for modelling molecular outflows. For example, the mass 
of H 2 outflowing with velocities > 50kms -1 is enhanced by a factor ~ 20 in non¬ 
equilibrium. We also compute the observable line emission from Cll and CO. Both 
are stronger at higher metallicity, while Cll and CO emission are higher for stronger 
and weaker radiation fields respectively. We find that Cll is generally unaffected by 
non-equilibrium chemistry. However, emission from CO varies by a factor of ^ 2 — 4. 
This has implications for the mean Xqo conversion factor between CO emission and 
H 2 column density, which we find is lowered by up to a factor ~ 2.3 in non-equilibrium, 
and for the fraction of CO-dark molecular gas. 

Key words: astrochemistry - ISM: atoms - ISM: molecules - galaxies: evolution - 
galaxies: ISM. 


1 INTRODUCTION 

Hydrodynamic simulations of galaxy formation typically 
model gas cooling by tabulating the cooling rate as a func¬ 
tion of gas properties such as the density and temperature, 
under certain assumptions. For example, the simplest ap¬ 
proach, as used in some of the first cosmological hydrody¬ 
namic simulations (e.g. Katz et al. 1992), is to assume that 
the gas has primordial abundances and is in collisional ion¬ 
isation equilibrium (CIE). Sutherland V Dopita (1993) also 
included the effects of metal-line cooling, computing cooling 
curves in CIE for a range of metallicities. 

Another effect that can be important for gas cooling is 
the presence of a photoionising UV radiation field, which can 
change the ionisation balance and heat the gas. Efstathiou 
(1992) showed that an extragalactic UV background (UVB) 
can suppress the cooling rate in a primordial plasma, thereby 
inhibiting the formation of dwarf galaxies. Katz et al. (1996) 


implemented primordial radiative cooling in the presence of 
a UVB in cosmological hydrodynamic simulations. 

Wiersma et al. (2009) considered the impact that a 
photoionising UVB has on cooling rates in the presence of 
metals. They showed that photoionisation can suppress the 
cooling rate by up to an order of magnitude at tempera¬ 
tures and densities typical of the intergalactic medium (e.g. 
10 4 K < T < 10 6 K, p/(p) < 100). They also showed that 
variations in relative abundances from their solar values can 
change the cooling rate by a factor of a few. Wiersma et 
al. (2009) tabulated the cooling rate from 11 elements sep¬ 
arately in the presence of the redshift-dependent UVB of 
Haardt V Madau (2001), and these tables have been used in 
several cosmological hydrodynamic simulations (e.g. Crain 
et al. 2009; Schaye et al. 2010, 2015; Hopkins et al. 2014). 

The effects of metal cooling and UV radiation are par¬ 
ticularly important below 10 4 K, as cooling from atomic hy¬ 
drogen becomes inefficient at such temperatures. In primor- 
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dial gas, the only major coolants are H 2 and HD (e.g. Saslaw 

V Zipoy 1967; Flower et al. 2000; Glover V Abel 2008). How¬ 
ever, these molecules can easily be dissociated by Lyman- 
Werner radiation (11.2eV — 13.6eV). If metals are present, 
efficient cooling can continue down to T < 100 K via fine- 
structure line emission (e.g. Glover V Jappsen 2007; Rich¬ 
ings et al. 2014a). Conversely, the presence of UV radiation 
will hinder cooling to such low temperatures, as it heats 
the gas via photoionisation and photoelectric heating from 
dust grains, the latter of which also depends on metallicity. 
Metallicity and UV radiation therefore influence the tran¬ 
sition from the warm (~ 10 4 K) to the cold 10 2 K) gas 
phase. 

Wolfire et al. (2003) investigated the warm-to-cold tran¬ 
sition by calculating the minimum pressure, and hence the 
minimum density, at which a cold phase can exist in pressure 
equilibrium with a warm phase, as a function of UV inten¬ 
sity, metallicity, dust abundance and ionisation rate from 
cosmic rays and Extreme UV (EUV)/X-ray radiation (see 
their equations 33-35). This minimum density increases with 
decreasing metallicity and increasing UV intensity (see also 
Glover & Clark 2014). 

Schaye (2004) demonstrated that the critical surface 
density above which a cold ISM phase can form in a galac¬ 
tic disc increases with decreasing metallicity and increasing 
UV intensity (see his equation 23). He additionally showed 
that the formation of a cold phase can trigger gravitational 
instabilities, due to the lower velocity dispersion in cold gas. 
Metallicity and UV radiation can thus affect the star form¬ 
ing properties of galaxies, in particular the threshold surface 
density below which star formation is unable to proceed. 
The transition from atomic to molecular hydrogen also de¬ 
pends on metallicity and UV radiation (e.g. Schaye 2001; 
Krumholz et al. 2008; Gnedin V Kravtsov 2011; Sternberg 
et al. 2014). 

Many of the above approaches assume that the gas is 
in chemical (including ionisation) equilibrium, or in other 
words, that the abundances of individual ions and molecules 
have reached an equilibrium or steady state. However, this 
assumption may not be valid if the dynamical or cooling 
time-scale of the gas is short compared to the chemical time- 
scale (e.g. Kafatos 1973; Sutherland & Dopita 1993; Gnat 
& Sternberg 2007; Oppenheimer & Schaye 2013a; Vasiliev 
2013), or if the UV radiation field is varying on time-scales 
shorter than the chemical time-scale (e.g. Oppenheimer V 
Schaye 2013b). 

Non-equilibrium chemistry has been shown to affect 
cooling rates in certain idealised scenarios. For example, in 
collisionally ionised gas that is cooling either isobarically 
or isochorically, the cooling rate computed self-consistently 
from non-equilibrium abundances is suppressed at temper¬ 
atures T > 10 4 K compared to gas in CIE (e.g. Sutherland 

V Dopita 1993; Gnat V Sternberg 2007). Oppenheimer V 
Schaye (2013a) investigated non-equilibrium ionisation and 
cooling in the presence of a photoionising radiation field, 
and they showed that the UVB reduces the non-equilibrium 
effects, although the cooling rates at T > 10 4 K are still 
suppressed with respect to chemical equilibrium. Richings 
et al. (2014a) showed that, at temperatures below 10 4 K, 
the direction of the non-equilibrium effects is reversed and 
the cooling rates are enhanced, due to an increase in the 
electron abundance compared to equilibrium. 


These studies considered special cases of a plasma that 
is cooling at constant density or pressure. Walch et al. 
(2011) compared equilibrium and non-equilibrium chem¬ 
istry and cooling in high-resolution simulations of a turbu¬ 
lent medium, in a box 500 pc across. They found that non¬ 
equilibrium chemistry did have a noticeable, albeit fairly 
small, effect on the gas temperature in their simulations. 
Vasiliev (2013) applied a treatment for non-equilibrium cool¬ 
ing to simulations of a supernova remnant, and showed that 
the evolution of the supernova remnant is sensitive to non¬ 
equilibrium effects. However, it remains to be seen whether 
such non-equilibrium effects will be relevant on galactic 
scales. 

In this paper we investigate the effects of metallicity, 
UV radiation and non-equilibrium chemistry on galaxy for¬ 
mation. We run a series of hydrodynamic simulations of iso¬ 
lated galaxies with stellar and total halo masses of 10 9 and 
10 11 M 0 , respectively. The simulations use a resolution of 
750 M© per particle and a gravitational force softening of 
3.1 pc, which is sufficient to resolve the warm-to-cold transi¬ 
tion. We consider a range of metallicities, 0.01 Z 0 ^ Z ^ Z 0 , 
and a range of UV radiation fields that span nearly three 
orders of magnitude in Hi photoionisation rate: the local 
interstellar radiation field (ISRF) of Black (1987), ten per 
cent of the Black (1987) ISRF and the redshift zero extra- 
galactic UV background (UVB) of Haardt & Madau (2001). 
We include a local prescription for self-shielding by gas and 
dust, and we also consider a model without self-shielding for 
comparison. We first run these simulations with the tem¬ 
perature and chemical abundances of the gas evolved us¬ 
ing the non-equilibrium chemical model of Richings et al. 
(20 14a, b), which follows the ionisation states of 11 elements 
that are important for the gas cooling rate, along with the 
abundances of 20 molecular species, including H 2 and CO. 
We then compare these to simulations evolved using cooling 
rates tabulated in chemical equilibrium. 

In addition to the dynamical impact on the galaxy due 
to the effects on the cooling rate, non-equilibrium chemistry 
can also affect observable diagnostics of individual chemi¬ 
cal species, e.g. in the presence of a fluctuating UV field 
(Oppenheimer & Schaye 2013b), or in the presence of su¬ 
personic turbulence (Gray et al. 2015). To investigate how 
non-equilibrium chemistry can affect observable emission on 
galactic scales, we perform radiative transfer calculations in 
post-processing to compute line emission from Cn and CO. 
We then compare the intensity of line emission computed 
from non-equilibrium abundances to that computed assum¬ 
ing chemical equilibrium. 

The remainder of this paper is organised as follows. In 
section 2 we describe the hydrodynamic methods and sub¬ 
grid physics models used in our simulations. We present our 
simulations in section 3, where we discuss the initial con¬ 
ditions (3.1), morphologies and star formation rates (3.2), 
outflow properties (3.3), and the phase structure of the ISM 
(3.4). We compute the observable line emission from Cn and 
CO in section 4, and we summarise our main results in sec¬ 
tion 5. 


MNRAS 000 , 1-26 (2016) 


Metallicity, radiation and chemistry in galaxies 3 


2 SUBGRID PHYSICS MODELS 

Our simulations were run using a modified version of the 
tree/Smoothed Particle Hydrodynamics (SPH) code GAD¬ 
GETS (last described in Springel 2005). The hydrodynamic 
equations were solved using the set of numerical meth¬ 
ods known as anarchy, which includes many of the lat¬ 
est improvements to the standard SPH implementation. In 
particular, anarchy uses the pressure-entropy formulation 
of SPH, derived by Hopkins (2013); the artificial viscosity 
switch of Cullen V Dehnen (2010); a switch for artificial 
conduction, similar to that used in Price (2008); the time- 
step limiters of Durier V Dalla Vecchia (2012); and the C 1 2 
Wendland (1995) kernel, anarchy will be described in more 
detail in Dalla Vecchia (in preparation); see also Appendix 
A of Schaye et al. (2015) for a description of the implementa¬ 
tion of anarchy in the cosmological simulations that were 
run for the eagle project, which is identical to the anarchy 
implementation that we use here. 

We use prescriptions to model physical processes that 
are unresolved in our simulations, including the chemical 
evolution of ions and molecules, radiative cooling, star for¬ 
mation and stellar feedback. These subgrid models are sum¬ 
marised below. 

2.1 Chemistry and radiative cooling 

We use the chemical model of Richings et al. (2014a,b) to 
evolve the chemical abundances of 157 species, including all 
ionisation states of the 11 elements that are most important 
for cooling 1 and 20 molecular species 2 . The gas temperature 
is evolved using radiative cooling and heating rates calcu¬ 
lated from the non-equilibrium abundances. This gives us 
a set of 158 differential equations (157 chemical rate equa¬ 
tions and the energy equation), which we integrate over each 
hydrodynamic timestep for each gas particle. We integrate 
these differential equations using the backward difference 
formula method and Newton iteration in cvode (from the 
sundials 3 suite of non-linear differential/algebraic equation 
solvers), using a relative tolerance of 10 -4 and an absolute 
tolerance of 10 -10 . 

We summarise the main chemical and thermal processes 
that are included in our model below. 


Table 1 . Properties of the UV radiation fields considered in this 
paper. 


UV field 

G 0 a 

r H , (s-U 

r H e, (s-T 

r H e„ (s-Y 

ISRF 

1.2 

4.4 x 10“ n 

3.7 x 10“ 12 

1.7 x 10“ 14 

lowISRF 

0.12 

4.4 x 10 -12 

3.7 x 10 -13 

1.7 x 10“ 15 

UVB 

0.014 

8.4 x 10“ 14 

2.0 x 10“ 14 

1.5 x 10“ 16 


a Radiation field strength in the energy range 6.0 — 13.6 eV, in 
units of the Habing (1968) field. 

^ Unattenuated Hi photoionisation rate. 
c Unattenuated Hei photoionisation rate. 

^ Unattenuated Hen photoionisation rate. 


(McElroy et al. 2013) where available, or using the equations 
from Lotz (1967), Silk (1970) and Langer (1978) otherwise. 
We include secondary ionisations of Hi and Hei from cosmic 
rays using the tables of Furlanetto & Stoever (2010), assum¬ 
ing a typical mean primary electron energy of E — 35 eV. 

We include photoionisation of atoms and ions, including 
Auger ionisation, using optically thin cross sections calcu¬ 
lated in the grey approximation for a given UV spectrum. 
We consider three UV spectra in this paper: the local ISRF 
of Black (1987), ten per cent of the Black (1987) ISRF and 
the redshift zero extragalactic UVB of Haardt V Madau 
(2001). The properties of these radiation fields are sum¬ 
marised in Table 1. We then attenuate the optically thin 
rates by the gas and dust, as a function of the column den¬ 
sities of Hi, H 2 , Hei, Hen and dust, as described in Richings 
et al. (2014b). To calculate these column densities, we as¬ 
sume that each gas particle is shielded locally. The total 
hydrogen column density, Vn tot , is then: 

N a tot= n a tot L, ( 2 . 1 ) 

where nH tot is the hydrogen number density of the gas par¬ 
ticle and L is the shielding length, i.e. the length scale over 
which the gas particle is able to shield itself. We use a 
Sobolev-like approximation to estimate the shielding length 
based on the gradient of the density, p : 


2.1.1 Chemical processes 

We include collisional ionisation, radiative and di-electronic 
recombinations and charge transfer reactions. The forma¬ 
tion of molecular hydrogen occurs on dust grains, using 
equation 18 of Cazaux V Tielens (2002) with a dust tem¬ 
perature Td U st = 10 K, as well as via gas phase reactions. 
We also include cosmic ray ionisations, assuming a primary 
ionisation rate of atomic hydrogen due to cosmic rays of 
Chi = 2.5 x 10 _17 s _1 (Williams et al. 1998). The primary 
ionisation rates of other species due to cosmic rays are then 
scaled to this value using the ratios in the umist database 4 


L-L Sob , p - | 2 £ p |. 

Gnedin et al. (2009) use this approximation in cosmological 
simulations to follow the formation of molecular hydrogen. 
They show that this approximation accurately reproduces 
the true column densities in their simulations, as measured 
along random lines of sight, with a scatter of a factor of 2 
in the range 3 X 10 20 cm -2 < Nhi + 2Nh 2 < 3 x 10 23 cm -2 . 
The column density of species i is then: 


1 H, He, C, N, O, Ne, Mg, Si, S, Ca & Fe. 

2 H 2 , H+, H+, OH, H 2 0, C 2 , 0 2 , HCO+, CH, CH 2 , CH+, CO ; 
CH+, CH+, OH+, H 2 0+, H 3 O+, CO+, HOC+, O+. 

3 https://computation.llnl.gov/casc/sundials/main.html 

4 http://www.udfa.net/ 


Ni — XiNu tot , ( 2 - 3 ) 

where Xi — m/rm tot is the abundance of species i. 

The photodissociation rates of molecular species are 
also attenuated using the shielding length given in equa- 
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tion 2.2. The attenuation of the photodissociation rate of 
species i due to dust is given by: 

Sd(Nn tot ,Z) = exp(—7dj4 v ), (2.4) 

where S^(Nu tot , Z) is the shielding factor (i.e. the ratio of 
the optically thick to optically thin photodissociation rates) 
due to dust, A v — 4.0 X 10 -22 Nn tot Z /Zq mag is the visual 
dust extinction, Z is the metallicity and the factors are 
taken from table 2 of van Dishoeck et al. (2006) where avail¬ 
able, or from table B2 of Glover et al. (2010) otherwise. 

In addition to dust shielding, molecular hydrogen can 
also be self-shielded. We use the temperature-dependent self¬ 
shielding function given in equations 3.12 to 3.15 of Richings 
et al. (2014b), which gives the shielding factor of H 2 as a 
function of the molecular hydrogen column density, 1Vh 2 , gas 
temperature, T, and Doppler broadening parameter, b. The 
Doppler broadening parameter includes thermal broadening 
and broadening due to turbulence. In our simulations, we 
assume a constant turbulent broadening parameter 5 tU rb — 

7.1 kms -1 , as used by e.g. Krumholz (2012), as this is typical 
for observed molecular clouds. 

CO can also be self-shielded and shielded by H 2 . We use 
the tables of the shielding factor of CO as a function of CO 
and H 2 column densities given by Visser et al. (2009). 

A full list of the chemical reactions that are included 
in our model can be found in table B1 of Richings et al. 
(2014a). 


2.1.2 Thermal processes 

We use the non-equilibrium abundances from the chemical 
network to calculate the net cooling rate of the gas. We in¬ 
clude cooling from the collisional excitation and ionisation of 
H and He, fine-structure line emission from metals, recombi¬ 
nation cooling and bremsstrahlung radiation. We use the H 2 
cooling function from Glover & Abel (2008) for rovibrational 
cooling from H 2 , and we include cooling from CO, H 2 O and 
OH. The photoheating rate is attenuated by gas and dust, as 
described by Richings et al. (2014b), using column densities 
calculated from equations 2.1 and 2.2 above. We also in¬ 
clude photoelectric heating on dust grains (Bakes & Tielens 
1994; Wolfire et al. 1995), cosmic ray heating (Goldsmith 
& Langer 1978; Glover & Jappsen 2007), and heating from 
the photodissociation of H 2 (Black & Dalgarno 1977), UV 
pumping of H 2 (Burton et al. 1990), and the formation of 
H 2 on dust grains (Hollenbach & McKee 1979) and in the 
gas phase (Karpas et al. 1979; Launay et al. 1991). 


2.1.3 Equilibrium cooling tables 

One aim of this paper is to compare simulations run using 
the full non-equilibrium chemical model described above to 
simulations evolved with cooling rates that are calculated 
assuming chemical equilibrium. We therefore used our chem¬ 
ical model to create tables of the cooling and heating rates 
and the mean molecular weight in chemical equilibrium as a 
function of gas temperature, hydrogen number density and 
total hydrogen column density. Note that the simulations in 
this paper were run at fixed metallicity, so we do not need to 


include an additional dimension for metallicity in the cool¬ 
ing tables. These tables were then used to evolve the gas 
temperature in the equilibrium cooling runs. By computing 
equilibrium cooling rates using the same chemical model as 
used for the non-equilibrium runs, we ensure that any dif¬ 
ferences are due to non-equilibrium effects, and not simply 
due to the use of different chemical rates in the models. 

2.2 Star formation 

We allow a gas particle to form stars if its temperature, T, 
and hydrogen number density, nH to t ? satisfy the following 
criteria: 

T < Tthresh = 1000 K, (2.5) 

riHtot > ^Htot,thresh = 1.0cm“ 3 . (2.6) 

Gas particles that satisfy these criteria form stars at a 
rate given by: 

p* = e (2.7) 

tff 

where p* is the star formation rate per unit volume, p gas 
is the gas density and ts = ^/37t/(32 Gp) is the free fall 
time. We use a star formation efficiency per free fall time 
of csf = 0.005. This is slightly lower than the value of 
csf ~ 0.015 that is typically observed in the Milky Way 
and nearby galaxies (Krumholz et al. 2012). We calibrated 
this parameter, along with the density and temperature 
thresholds and the parameters for the stellar feedback model 
described in the next section, to reproduce the observed 
Kennicutt-Schmidt relation in nearby galaxies (see Fig. 3). 

In a timestep At, star forming gas particles are stochas¬ 
tically turned into star particles with a probability p given 
by: 

p = min(^t l) . (2.8) 

V Pgas J 

2.3 Stellar feedback 

Each star particle represents a stellar population, rather 
than an individual star. We assume that the stellar pop¬ 
ulation initially follows the Chabrier (2003) Initial Mass 
Function (IMF) with masses in the range 0.1 — 100 M 0 . We 
then calculate how many type II supernovae will explode in 
each timestep for each star particle, using the metallicity- 
dependent stellar lifetimes of Portinari et al. (1998) and as¬ 
suming that all stars with a mass greater than 6M 0 will 
end their lives in a supernova (stars with masses of 6 — 8 M 0 
explode as electron capture supernovae in models with con¬ 
vective overshoot, e.g. Chiosi et al. 1992). 

The energetic feedback from supernovae is implemented 
using the stochastic thermal feedback prescription of Dalla 
Vecchia & Schaye (2012), except that we distribute the to¬ 
tal available energy from type II supernovae from a single 
star particle in time according to the stellar lifetimes of the 
massive stars, rather than combine it into a single supernova 
event. For each star particle that has a non-zero number of 
supernovae, TVsnii, in timestep At, we stochastically select 
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gas particles from the 7V ng b = 48 neighbours to be heated 
by AT = 10 7 ' 5 K. Note that N ng b is smaller than the num¬ 
ber of neighbours that we use to compute the SPH kernel 
(^nTb H = 100). This reduces the computational cost of the 
neighbour-finding routine for star particles, and has little 
impact on the accuracy of the stellar feedback. The proba¬ 
bility p of selecting a given gas particle to be heated is: 


p — min 


2 EsmiNsmipm p 

3 fc B ATEtT b mi ’ 



(2.9) 


where Esmt = 10 51 erg is the energy injected per supernova, 
p is the mean molecular mass of the gas particle, m p is the 
mass of a proton and rm is the mass of the i th neighbouring 
gas particle. 

By imposing a minimum heating temperature AT, we 
ensure that the cooling time of the heated gas is long enough 
to avoid artificial radiative losses that would otherwise make 
the feedback scheme ineffective. However, note that, at the 
resolution that we use in our simulations (750 Mq per par¬ 
ticle), we require approximately 10 supernovae to heat a 
single gas particle to 10 7 ' 5 K. We therefore do not resolve 
individual supernovae, as each feedback event is equivalent 
to several supernovae exploding simultaneously. 

While we include the energetic feedback from super¬ 
novae, we do not include the chemical enrichment from su¬ 
pernovae or from stellar mass loss. This allows us to follow 
model galaxies at a fixed metallicity, which eases the inter¬ 
pretation of our numerical experiments. 


2.4 Jeans limiter 

To ensure that we always resolve the Jeans mass, Mj, 
and the Jeans length, Lj, we impose a minimum floor 
on the pressure that enters the hydrodynamic equations. 
This is similar to the approach used by e.g. Robertson & 
Kravtsov (2008); Schaye V Dalla Vecchia (2008); Hopkins 
et al. ( 2011 ). 

There are two criteria that we can consider to determine 
the pressure floor. Firstly, we can require that the Jeans 
mass be resolved by at least a factor Nj jin times the mass 
within the SPH kernel. The Jeans mass is given by: 


Mj 


AAA 

GG^p 1 / 2 ’ 


( 2 . 10 ) 


where c s — is the sound speed, 7 = 5/3 is the ratio 

of specific heats, P is the thermal gas pressure and p is the 
gas density. 

If each gas particle has mass ra gas , and we use N^ g \ b H 
neighbours in the SPH kernel, then the mass within the 
kernel is: 


M k - Nngb m g as• (2.11) 

From equations 2.10 and 2.11, we find that the minimum 
pressure that we require is: 

Pfloor.m = (5) ' ^(iVj,miVn S g P b H "lgas) 2/ V /3 - (2-12) 


Secondly, we can require that the Jeans length be re¬ 
solved by at least a factor 7Vj,i times the SPH smoothing 
length, hsmi• Using equation 2.10, the Jeans length is: 


T _ s 

J - G 1 / 2 /? 1 / 2 ' 

Hence the minimum pressure that we require is: 


(2.13) 


Pfloor,i= — JV|,ifcLiP a . (2-14) 

7ry 

Since the smoothing length is defined such that the 
number of SPH neighbours, N .is constant, we can ex¬ 
press h sm \ in terms of the kernel mass as: 


hsmi — 


3 iv; 


ngb 


47rp 


1/3 


(2.15) 


By combining equations 2.12, 2.14 and 2.15, we see that 
these two criteria are equivalent if: 


Vm -2 (2-16) 

We set the pressure floor in our simulations such that 
we resolve Mj by at least Nj : m = 4 kernel masses, and thus 
we resolve Tj by at least IVj,i = 3.2 smoothing lengths. 

Note that, while we impose this floor on the pressure 
that enters the hydrodynamic equations, the temperature 
is still allowed to cool below this limit. This means that 
the thermal pressure will effectively be decoupled from the 
hydrodynamic equations once this floor is reached. For the 
resolution of our simulations, the pressure floor that we use 
corresponds to the thermal pressure at a temperature of 
240 K at the star formation threshold density nH tot ,thresh = 
1.0 cm -3 . We implement the Jeans limiter as a pressure floor 
rather than a temperature floor (as was used in Schaye & 
Dalla Vecchia 2008) so that the thermal and chemical state 
of the gas will evolve towards a realistic equilibrium for 
the given density, although we miss small-scale, high-density 
structures that are unresolved in the simulation. 


3 SIMULATIONS 

We ran a suite of SPH simulations of isolated galaxies with 
a range of metallicities and UV radiation fields, using a res¬ 
olution of 750 M® per gas particle, with 100 SPH neigh¬ 
bours, and a gravitational softening length of 3.1 pc. Each 
simulation was evolved for 1 Gyr. We ran each simulation 
twice: once using the full non-equilibrium chemical model 
summarised in section 2 . 1 , and once using tabulated cooling 
rates calculated assuming chemical equilibrium, as described 
in section 2.1.3. In this section, we describe the initial con¬ 
ditions, and we compare the morphologies, star formation 
rates, outflow properties and ISM phases in our simulated 
galaxies. For each of these, we focus on the effects of metal¬ 
licity, radiation field and non-equilibrium chemistry. 
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Table 2. Parameters and properties of the galaxies in the suite of simulations used in this paper: total mass M 200 within the 
radius R200,crit enclosing a mean density of 200 times the critical density of the Universe at redshift zero, concentration C 200 of 
the dark matter halo, initial stellar mass initial gas mass M gaS; i n i t , disc gas mass fraction /d, ga s? mass per gas or star 

particle mbaryon? gravitational softening length e so f t , gas metallicity Z, UV radiation field and whether we include self-shielding. 


Model 

M 200 

(M 0 ) 

C200 

, init 

(M 0 ) 

A^gas, init 

(M 0 ) 

/d,gas 

^baryon 

(M®) 

e soft 

(pc) 

z 
(Z ©) 

UV Field a 

Shielding 

ref 

10 11 

8.0 

1.4 x 10 9 

4.8 x 10 s 

0.3 

750 

3.1 

0.1 

ISRF 

yes 

lowZ 

10 11 

8.0 

1.4 x 10 9 

4.8 x 10 8 

0.3 

750 

3.1 

0.01 

ISRF 

yes 

hiZ 

10 11 

8.0 

1.4 x 10 9 

4.8 x 10 s 

0.3 

750 

3.1 

1.0 

ISRF 

yes 

lowISRF 

10 11 

8.0 

1.4 x 10 9 

4.8 x 10 8 

0.3 

750 

3.1 

0.1 

lowISRF 

yes 

UVB 

10 11 

8.0 

1.4 x 10 9 

4.8 x 10 s 

0.3 

750 

3.1 

0.1 

UVB 

yes 

UVBthin 

10 11 

8.0 

1.4 x 10 9 

4.8 x 10 8 

0.3 

750 

3.1 

0.1 

UVB 

no 


a See Table 1 

3.1 Initial conditions 

The initial conditions that we use are based on the model of 
Springel et al. (2005), and were generated using a modified 
version of a code that was kindly provided to us by Volker 
Springel. Each galaxy consists of an exponential disc of gas 
and stars with a radial scale length of 2.0 kpc, and a central 
stellar bulge with a Hernquist (1990) density profile, embed¬ 
ded within a dark matter halo. The stellar disc and bulge are 
represented by collisionless particles with the same mass as 
the gas particles (750M© per particle). The main difference 
between our initial conditions and the model of Springel et 
al. (2005) is that we represent the dark matter halo using a 
static potential, rather than with live dark matter particles. 
Using a static potential speeds up the calculation without 
affecting the results. 

The gas is initially isothermal, with a temperature of 
10 4 K, and the chemical abundances are initially in chem¬ 
ical equilibrium. The stellar disc is set up with a vertical 
distribution that follows an isothermal profile with a scale 
height of ten per cent of the radial scale length of the disc, 
while the vertical structure of the gaseous disc is set to be 
in hydrostatic equilibrium using an iterative procedure. 

The total mass of each galaxy within i?2oo,crit (i.e. 
the radius enclosing a sphere with a mean density of 200 
times the critical density of the Universe at redshift zero) 
is M200 = 1O 11 M 0 , and the initial stellar mass is M* = 
1.4 x 1O 9 M 0 . These masses are consistent with the stellar 
mass-halo mass relation obtained from abundance matching 
by Moster et al. (2013) and corrected for baryonic effects 
using the prescription of Sawala et al. (2015). 

A fraction /* s b = 0.2 of the stellar mass is in the bulge, 
with the remainder in the disc. We use a disc gas mass 
fraction /d, gas = 0.3, which gives an initial gas mass of 
M gas = 1.8 x 10 8 M 0 . The dark matter density profile that 
we use to calculate the static dark matter potential follows 
a Hernquist (1990) profile that has been scaled to match the 
inner regions of a Navarro et al. (1996) (NFW) profile with 
a concentration C 200 = 8.0, which agrees with the redshift 
zero mass-concentration relation of Duffy et al. (2008). 

In our reference model (ref), we use a fixed metallicity 
of 0.1 Z 0 , and for the UV radiation field we use the ISRF 
of Black (1987), along with the self-shielding prescription 
described in section 2.1.1. We also consider two additional 
metallicities (0.01 Z 0 and Z 0 ), and three additional radia¬ 
tion fields (ten per cent of the Black (1987) ISRF and the 


redshift zero UVB of Haardt V Madau (2001), both with 
self-shielding, and the redshift zero Haardt & Madau (2001) 
UVB without self-shielding). The parameters and properties 
of the galaxies in our suite of simulations are summarised 
in Table 2. We summarise the properties of the different 
radiation fields used in our simulations in Table 1. 

Throughout this paper we assume that the dust-to-gas 
ratio scales linearly with metallicity, or in other words, that 
the dust-to-metals ratio is constant. However, there is ob¬ 
servational evidence that the dust-to-metals ratio decreases 
at low metallicity, below ~ O.3Z 0 (e.g. Remy-Ruyer et al. 
2014). Therefore, it is possible that our runs at 0.01 Z 0 and 
0.1 Z Q overestimate the total dust abundance. This would 
affect the formation rate of H 2 , shielding of the radiation 
field, and photoelectric heating from dust grains. 

3.2 Morphology and star formation 

Maps of the gas surface density after 500 Myr from the sim¬ 
ulations run with the full non-equilibrium chemical model 
are shown in Fig. 1 for different metallicities (top row) and 
different UV radiation fields (bottom row), and maps of the 
gas temperature after 500 Myr are shown in Fig. 2. Each 
pair of panels in these figures shows projections looking at 
the disc face-on (top) and edge-on (bottom). 

Comparing the three different metallicities, we see that 
the morphology of the gas is very different in these three 
runs. In the simulation with the lowest metallicity (one per 
cent solar; top left pair of panels in Figs. 1 and 2), star 
formation only occurs at the centre of the disc. However, at 
higher metallicities, star formation becomes more vigorous 
and extends to larger radii. This leads to more supernovae, 
which create more bubbles of hot gas in the disc (as seen 
in Fig. 2), and drive more gas out of the disc in vertical 
fountains and outflows (as seen in the edge-on views). 

These trends with metallicity occur because, at higher 
metallicities, there is more metal-line cooling, which allows 
the gas to cool down to a cold ISM phase (with temperatures 
of a few hundred Kelvin) at lower densities. Gas needs to 
cool down to the cold ISM phase before it can form stars. 
Therefore, in the runs at higher metallicity, star formation 
can proceed in lower-density gas, and hence at lower gas 
surface densities. Our models all start with the same gas 
density profile. Hence, if star formation extends to lower gas 
surface densities, then it will extend to larger radii. 
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Figure 1. Maps of the total gas surface density after 500 Myr in the simulations run with the full non-equilibrium chemical model. The 
simulations in the top row were run with the Black (1987) ISRF at different metallicities: 0.01 Zq (lowZ; left), 0.1 Z© (ref; centre) and 
Z© (hiZ; right). The simulations in the bottom row were run at a metallicity of 0.1 Z© with different UV radiation fields: ten per cent 
of the Black (1987) ISRF (lowISRF; left), the Haardt & Madau (2001) UVB at redshift zero ( UVB; centre) and the Haardt Sz Madau 
(2001) UVB without self-shielding (UVBthin; right). Each pair of panels shows projections looking at the disc face-on (top) and edge-on 
(bottom). The face-on projections are 8kpc across, and the edge-on projections cover 4kpc in the vertical direction. 


Additionally, the increase in metal-line cooling at higher 
metallicities means that gas can more easily undergo gravita¬ 
tional collapse to higher densities. These two effects (star for¬ 
mation at lower densities and more gas collapsing to higher 
densities) lead to an increase in the total star formation rate, 
and hence more stellar feedback in the disc. 

Furthermore, the fact that the Jeans scale is smaller at 
lower temperatures also contributes to the more fragmented 
appearance of the gas disc at high metallicities. 

However, it is important to note that these runs at dif¬ 
ferent metallicities all use the same radiation field. In reality, 
the lower star formation rate that we find at lower metallic¬ 
ity will also result in a weaker radiation field. As we show 
below, this will tend to increase the star formation rate and 
thus will lessen the differences in star formation rate at dif¬ 
ferent metallicities. 

The bottom rows of Figs. 1 and 2 show runs at ten per 
cent solar metallicity with different UV radiation fields: ten 
per cent of the Black (1987) ISRF (‘lowISRF’) and the red- 
shift zero extragalactic UVB of Haardt & Madau (2001) 
(‘UVB’), both run with self-shielding, and the redshift 
zero Haardt Madau (2001) UVB without self-shielding 


(‘UVBthin’). See Table 1 for the properties of these radiation 
fields. As we decrease the strength of the UV radiation field 
(from the reference run in the top centre panels of Figs. 1 
and 2, to ‘lowISRF’ to ‘UVB’), the gas disc becomes more 
fragmented. A lower radiation field reduces the heating rate 
from photoionisation and photoelectric dust heating, which 
allows the gas to cool to the cold, star forming ISM phase 
(with temperatures of a few hundred Kelvin) at lower densi¬ 
ties. This increases the star formation rate, leading to more 
stellar feedback, and it decreases the Jeans scale in this gas. 
Hence the gas becomes more fragmented, similar to the ef¬ 
fect of increasing the metallicity. 

In the run without self-shielding (‘UVBthin’; bottom 
right pair of panels in Figs. 1 and 2), the gas is heated 
by photoionisation of Hi even at high densities, where it 
would otherwise become shielded from ionising radiation. 
This means that in this run gas can only cool below 1000 K 
and form stars at high densities (riH,tot ^ 10cm -3 ). Fur¬ 
thermore, the high-density gas remains much warmer than 
in the corresponding run with self-shielding (‘UVB’), so the 
Jeans scale in this gas is larger. Hence the gas morphology 
is smoother when we do not include self-shielding. 
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Figure 2. As Fig. 1, but for the gas temperature. Except for Z = 0.01 Zq, hot bubbles of gas driven by supernovae disrupt the otherwise 
smooth density distribution in the disc and drive vertical outflows of gas. 


The simulations in Figs. 1 and 2 were run with the full 
non-equilibrium chemical model. In the simulations evolved 
with cooling rates in chemical equilibrium, we find very sim¬ 
ilar morphologies to those seen in Figs. 1 and 2. 

The trends of star formation rate with metallicity and 
radiation field can be understood more clearly in plots of 
the star formation rate surface density, Esfr, versus gas 
surface density of neutral (atomic and molecular) hydro¬ 
gen, Ehi+h 2 7 i-e. the Kennicutt-Schmidt relation (Kennicutt 
1998). These are shown in Fig. 3 for different metallicities 
(top row) and different radiation fields (bottom row). The in¬ 
tensity of the blue-green colour indicates the mass-weighted 
distribution of gas as a function of Esfr and Ehi+h 2 in 
simulations evolved with the full non-equilibrium chemical 
model. We measure these on a two-dimensional grid of cells, 
each 100 pc across, that span a square region 8 kpc across 
centred on the disc and are aligned such that the disc is 
viewed face on. We include all particles within ±2 kpc of the 
mid-plane of the disc in the vertical direction. We calculate 
Esfr from the gas particles in the cell. For each gas parti¬ 
cle that meets our star formation criteria (nH tot > 1 cm -3 
and T < 1000 K), we compute the star formation rate using 
equation 2.7 for the star formation rate density, multiplied 
by the volume of the particle (m gas /pgas)- We then sum the 
star formation rates of all star-forming particles in the cell 


to obtain the total instantaneous star formation rate, which 
we divide by the area of the cell to get Esfr. 

In each simulation we have combined measurements of 
Esfr and Ehi+h 2 from all snapshot outputs from 100 Myr 
to 1000 Myr, at intervals of 100 Myr. The cyan curves in 
Fig. 3 show the median value of Esfr in bins of Ehi+h 2 
for simulations evolved with the full non-equilibrium model 
(‘NonEq’; solid curves) and with cooling rates in chemical 
equilibrium (‘Eqm’; dashed curves), and the contours show 
the observed Kennicutt-Schmidt relations in various samples 
of galaxies, as indicated in the legend. 

At high gas surface densities, Esfr tends towards a 
power-law relation, with a slope similar to (albeit slightly 
larger than) that measured by Kennicutt (1998) (the red 
region in Fig. 3). At low gas surface densities, the star for¬ 
mation rate drops below this power-law relation, and is cut 
off below a threshold gas surface density, Ehi+h 2 • 

In the top right panel of Fig. 3, we see that the simula¬ 
tion at solar metallicity agrees fairly well with the relation 
measured for nearby galaxies by Bigiel et al. (2008) (black 
and grey contours). Note that we calibrated the parameters 
of our star formation and stellar feedback models to match 
the normalisation of the observed relation. The threshold 
E“ below which star formation is cut off increases with 
decreasing metallicity, as predicted by Schaye (2004). If we 
define Ehi+h 2 as the surface density at which the median 
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Figure 3. The star formation rate surface density, £sfr> versus surface density of neutral (atomic plus molecular) hydrogen, Ehi+h 2 j be. 
the Kennicutt-Schmidt relation. The intensity of the blue-green colour corresponds to two-dimensional histograms of Esfr and Ehi+h 2 
measured in the simulations run with the full non-equilibrium chemical model at different metallicities (top row) and for different radiation 
fields (bottom row; see Table 1). We measure Esfr and Ehi+h 2 on a grid of cells, each 100 pc across, that span a square region 8kpc 
across aligned such that the disc is viewed face on. We combine measurements from all snapshot outputs from 100 Myr to 1000 Myr, 
at intervals of 100 Myr. We set all cells to have a minimum E$fr of 10 -6 M@ yr _1 kpc -2 , so that they are visible in these plots. The 
cyan curves show the median value of E$fr in bins of Ehi+h 2 , from simulations evolved with the full non-equilibrium chemical model 
(solid) and those run with cooling rates in chemical equilibrium (dashed). The contours show the observed Kennicutt-Schmidt relation 
from various samples. The red shaded region shows the best-fit power-law relation of Kennicutt (1998) (their equation 4), with the 
width of this region indicating the uncertainty in the normalisation. The black and grey contours are taken from Bigiel et al. (2008) (the 
centre-right panel of their fig. 8), measured in sub-kpc regions of nearby galaxies. Finally, the yellow contours are take from Bolatto et 
al. (2011) (from their fig. 6), measured in the Small Magellanic Cloud. We use their data at a resolution of 200pc, rather than their full 
resolution data, as this is closer to the spatial scale that we use for the measurements of Ssfr an d ^hi+h 2 in the simulations (100pc). 
We have renormalised SgpR from the observations for a Chabrier (2003) IMF, as used in our simulations. 
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Esfr is 1 dex below the power-law relation of Kennicutt 
(1998), we find ^m r +H 2 oc Z~ 03 . This agrees well with the 
metallicity dependence derived by Schaye (2004) for the col¬ 
umn density at which the cold ISM phase forms (see his 
equation 23). 

In the lowest-metallicity simulation, in the top left panel 
of Fig. 3, the gas surface densities do not extend much above 
E“ 5 so we cannot see whether the relation turns over and 
follows a power-law relation at higher gas surface densities. 

Comparing runs at ten per cent solar metallicity for dif¬ 
ferent radiation fields (top centre, bottom left and bottom 
centre panels of Fig. 3), we see that, as we decrease the 
strength of the UV radiation field, the threshold gas sur¬ 
face density, Ehi+h 2 , below which star formation is cut off 
decreases, as predicted by Schaye (2004). This is similar to 
the effect of increasing the metallicity that we see in the top 
row. We find Ehi+h 2 oc Gq ' 3 , in good agreement with the 
dependence on UV intensity derived by Schaye (2004) for 
the column density at which the cold ISM phase forms (his 
equation 23). 

The bottom right panel shows that, when we do not in¬ 
clude self-shielding, Ehi+h^ increases. Furthermore, at gas 
surface densities above Ehi+hv, , the star formation rate con¬ 
tinues to rise steeply, and lies above the observed relation of 
Kennicutt (1998). However, since we calibrated the param¬ 
eters of our star formation and stellar feedback models to 
reproduce the observed Kennicutt-Schmidt relation in sim¬ 
ulations that included self-shielding, this could possibly be 
remedied by re-calibrating these parameters. 

The star formation histories are shown in Fig. 4 for dif¬ 
ferent metallicities (top panel) and different radiation fields 
(bottom panel). The solid and dashed curves show simula¬ 
tions evolved with the full non-equilibrium chemical model 
and with cooling rates in chemical equilibrium respectively. 
The star formation rate initially rises rapidly as gas cools 
from its initial temperature of 10 4 K to form a cold, star¬ 
forming ISM phase. Once stellar feedback takes effect (af¬ 
ter ~ 50Myr), the star formation rate levels off and then 
steadily declines over the course of the simulation, as gas is 
either consumed by star formation or driven out of the disc 
in outflows. 

The star formation rate increases with metallicity and 
is about two orders of magnitude higher for solar metallicity 
than for 0.01 Z®. As we decrease the strength of the radia¬ 
tion field, the total star formation rate increases, by a factor 
of ~ 3 for the radiation fields that we consider here. This is 
similar to the trend that we see when we increase the metal¬ 
licity from 0.01 Z© to Z®, although the size of the effect is 
smaller than when we vary the metallicity over this range. 

In the run without self-shielding, the total star forma¬ 
tion rate is typically lower, by up to an order of magnitude, 
than in the corresponding run with self-shielding, although 
they are similar (to within a factor of two) between 400 and 
600 Myr. 

Comparing the solid and dashed curves in Figs. 3 and 
4, we see that non-equilibrium cooling has no noticeable 
systematic effect on the simulated Kennicutt-Schmidt re¬ 
lation or on the total star formation rate of the simulated 
galaxy. However, this conclusion may be dependent on the 
resolution of our simulations. In particular, our star for¬ 
mation prescription allows gas to form stars at densities 
tih > 1.0 cm -3 and temperatures T < 1000 K. In other 



Time (Myr) 


Figure 4. Star formation histories of simulations run with the 
full non-equilibrium chemical model (solid curves) and run with 
cooling rates in chemical equilibrium (dashed curves) at different 
metallicities (top panel) and for different radiation fields (bottom 
panel; see Table 1 ). The black curves in the two panels are from 
the same simulation. We find higher star formation rates in the 
simulations at higher metallicity and, to a lesser extent, in the 
presence of a weaker UV radiation field. 


words, gas in our models becomes star forming once it transi¬ 
tions from the Warm Neutral Medium (WNM) to the Cold 
Neutral Medium (CNM). It is possible that there are still 
important non-equilibrium effects on smaller scales than we 
resolve that could affect the star formation rate. 


3.3 Outflows 

Stellar feedback drives gas out of the disc in our simulations. 
To measure the radial mass outflow rates and velocities, we 
consider a spherical shell of radius 0. 2 ^ 200 ,crit = 19kpc and 
thickness Ar = i^2oo,crit/150 = 0.64 kpc, centred on the ori¬ 
gin (which is defined as the centre of the static dark mat¬ 
ter potential). The net mass outflow rate, dM net /dt , is then 
given by: 


dM net 

dt 


1 

A r 


jy-shell 



(3.1) 


where m*, v* and are the mass, velocity and position of 
the i th particle respectively, and we sum over the iV she11 gas 
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Figure 5. Evolution of the net mass outflow rate (left panels), mass loading factor (centre panels) and mass-weighted mean radial 
outflow velocity (right panels) measured through a thin spherical shell of radius 0.2i?200,crit — 19kpc and thickness 0.64 kpc. We show 
simulations evolved using the full non-equilibrium chemical model (solid curves) or using cooling rates in chemical equilibrium (dashed 
curves) at different metallicities (top row) and for different radiation fields (bottom row; see Table 1). The mass outflow rate increases 
with increasing metallicity and decreasing radiation field, although the mass loading factor at late times is ~ 10, regardless of metallicity 
or radiation field. 


particles in the shell. We also calculate the mass-weighted 
mean radial velocity, v ou t, of the outflowing particles as: 


Aj-shell 

_, x _£t + i m< ( v » • r */Jr* |) 

^out — W/+ — jyshell 

£i= + i m 


(3.2) 


where the subscript c +’ indicates that we only consider par¬ 
ticles with Vi • r i > 0, i.e. that are moving radially outward. 

We could split the net outflow rate in equation 3.1 into 
separate outflow and inflow rates, using particles that are 
moving radially outward or inward respectively. However, 
we find that the inflow rates at this radius are negligible in 
our simulations. 

We measure the mass outflow rates and mean outflow 
velocities from each simulation in snapshots output at 1 Myr 
intervals, using equations 3.1 and 3.2. These are shown in 
Fig. 5 for the simulations with different metallicities (top 
row), and in the presence of different UV radiation fields 
(bottom row). Solid and dashed curves are from simulations 
run with the full non-equilibrium chemical model and cool¬ 
ing rates in chemical equilibrium respectively. 

We see that the mass outflow rates (left panels) gener¬ 
ally increase during the first 200 Myr, as there is a delay 
between the onset of star formation and the first super¬ 
novae, and it takes a finite time for the gas to reach the 
radius of 19 kpc where we measure the outflows (~ 10 Myr 
for gas travelling at 200kms _1 ). In the simulations without 
self-shielding (green curves in the bottom row), this initial 
increase is more gentle and extends over a longer period of 
time (~ 500Myr). We saw in Fig. 4 that the star formation 
rates in the simulations without self-shielding also increase 
more gently over this period. After the initial rise, the mass 


outflow rates fluctuate around an approximately steady or 
gently declining value. 

The outflow rates tend to be larger in simulations at 
higher metallicity or in the presence of a weaker UV radi¬ 
ation field, due to the larger star formation rates that we 
find in these cases (see Fig. 4). In the centre panels of Fig. 5 
we show the evolution of the ratio between the mass outflow 
rate and the star formation rate, i.e. the mass loading fac¬ 
tor. After 200 Myr, the mass loading factor tends towards a 
value of rsj 10, independent of the metallicity or strength of 
the UV radiation field. An exception to this trend is seen 
in the simulations that do not include self-shielding (green 
curves in the bottom row), which show a steady rise in the 
mass loading factor from ~ 1 at 200 Myr to ~ 30 at the end 
of the simulation, albeit with large fluctuations on timescales 
of a few Myr. 

One caveat is that we do not include a gaseous halo in 
the initial conditions. In a realistic galaxy, the outflowing 
gas may be slowed down as it interacts with an existing 
gaseous halo, which could decrease the mass loading factor 
at i22oo,crit- Alternatively, it may sweep up more material 
from the halo, which could increase the mass loading factor. 

In the right panels of Fig. 5, we show the evolution of 
the mean radial velocity of outflowing gas. Initially, we find 
very large outflow velocities (v out ~ 1000 km s -1 ), because 
the fastest particles ejected from the disc are the first to 
reach the radius at which we measure the outflows. The 
mean outflow velocity subsequently decreases as the slower 
particles reach this radius, and tends towards a value of 
~ 65kms _1 by the end of the simulation. 

Comparing solid and dashed curves in Fig. 5, we see that 
the outflows are generally not strongly affected by the use 
of equilibrium cooling rates. In the runs at solar and ten per 
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cent solar metallicity in the presence of the Black (1987) 
ISRF (red and black curves respectively in the top row of 
Fig. 5), the mass outflow rates (left panels) rise more gently 
in the first 200 Myr in the simulations run with equilibrium 
cooling rates (dashed curves), compared to the correspond¬ 
ing runs evolved with non-equilibrium cooling rates (solid 
curves). However, this difference becomes less pronounced 
in the presence of a weaker UV radiation field (red and blue 
curves in the bottom row of Fig. 5). Also, the mass loading 
factors (centre panels) and outflow velocities (right panels) 
are similar whether we use non-equilibrium or equilibrium 
cooling. 

We also considered how the mass loading factor, (3 = 
M out /M*, depends on gas surface density in our simulations. 
To do this, we measured vertical outflows close to the disc, 
at a vertical height of 1 kpc above and below the mid-plane 
of the disc. We positioned two thin sheets of thickness A z = 
50 pc parallel to the disc (in the x — y plane) at a vertical 
distance \z\ = 1 kpc. Each sheet covered a square region 
8 kpc across centred on the galaxy centre, which corresponds 
to the region shown in the face-on views in Figs. 1 and 2. 
The vertical mass outflow rate through these two sheets is 
given by: 


dMo U t 

dt 


1 y- 

i=1 


(3-3) 
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Figure 6. Mass loading factor, (3 = M ou t/M*, plotted against 
gas surface density, E gas , measured in regions of the disc 100 pc 
across at a vertical distance \z\ of 1 kpc above and below the 
mid-plane of the disc. We show simulations at different metallic- 
ities (top row) and for different radiation fields (bottom row; see 
Table 1 ), evolved using the full non-equilibrium chemical model 
(black) or using cooling rates in chemical equilibrium (red). The 
solid curves show the median mass loading factor in bins of E gas , 
and the shaded regions indicate the range between the tenth and 
ninetieth percentiles in each bin. The mass loading factor de¬ 
creases with increasing E gas , and this relation is unaffected by 
using equilibrium cooling. 


where v z ,i is the ^-component of the velocity of the i th parti¬ 
cle, the subscript c +’ indicates that we only include particles 
with z • v > 0 (i.e. that are moving vertically away from the 
mid-plane), and we sum over the iV+ heet gas particles with 
z • v > 0 in the two sheets. 

To compare the mass loading factor at a vertical height 
of 1 kpc to the gas surface density within the disc, we divided 
the two thin sheets into a grid of cells, each 100 pc across. 
We then measured the mass outflow rates through the two 
sheets in each cell, along with the star formation rate and 
gas surface density in the disc within the cell. 

In Fig. 6 we plot the mass loading factor as a function of 
gas surface density for different metallicities (top row) and 
different UV radiation fields (bottom row). We measured the 
mass loading factors and gas surface densities in all snapshot 
outputs from each simulation, at intervals of 1 Myr, and we 
show the median mass loading factor (solid curves) and the 
range between the tenth and ninetieth percentiles (shaded 
regions) in bins of gas surface density. Simulations run with 
the full non-equilibrium and equilibrium chemical models 
are shown in black and red respectively. 

The mass loading factor decreases with increasing gas 
surface density, following approximately a power-law, with 
/3 oc E ga 2 s . At constant gas surface density, the mass loading 
factor increases with decreasing metallicity and increasing 
UV radiation field. The mass loading factors in the different 
panels of Fig. 6 suggest a scaling of (3 oc Z -0 ' 5 and (3 oc 
Gq' 3 , although these scalings are highly uncertain, as we 
only consider a small number of different metallicities and 
radiation fields. 

Comparing black and red curves in Fig. 6, we see that 
we recover the same relation between mass loading factor 
and gas surface density, regardless of whether we use non¬ 


equilibrium cooling rates or cooling rates in chemical equi¬ 
librium. 

Hopkins et al. (2012) explored the depen¬ 
dence of [3 on galaxy properties in their SPH 
simulations of isolated galaxies. They found = 

10(Vc(-R)/100kms -1 ) -1 (E g as(R)/10M G pc -2 ) -0 ' 5 (their 
equation 8), where Vc is the circular velocity at radius 
R. This has a much weaker dependence on E gas , with 
a power-law slope of —0.5, compared to ~ — 2 from our 
simulations. This different scaling with E gas may be due to 
the different prescription for stellar feedback that was used 
by Hopkins et al. (2012). 

The relation betweem [3 and E gas has also been explored 
by Creasey et al. (2013). They used high-resolution mesh 
simulations of supernova-driven galactic winds in a 1 kpc col¬ 
umn through a galactic disc, but they did not include radia¬ 
tive cooling below 10 4 K. They showed that /3 increases with 
decreasing E gas , and with increasing gas fraction, / gas . They 
fit a power-law dependence of (3 on E gas and / gas in their 
simulations, and they found (3 — 13E gas 15 / ga g 6 (see equa¬ 
tions 39-42 of Creasey et al. 2013). This power-law scaling 
with E gas is intermediate between our value and the relation 
of Hopkins et al. (2012). 

We find much larger mass loading factors than Creasey 
et al. (2013) at a given gas surface density. For example, at 
E gas = 10 Mq pc -2 , we typically find a median mass loading 
factor of ~ 100 at \z\ = lkpc, whereas, for a gas fraction 
/gas = 0.3, the best-fit relation of Creasey et al. (2013) gives 
a mass loading factor of 0.8. For comparison, the relation of 
Hopkins et al. (2012) gives /3 = 20, where Vc ~ 50kms _1 in 
our simulations. Creasey et al. (2013) measured the outflows 
from their simulation volume closer to the mid-plane, at 
\z\ = 500 pc. However, when we consider outflows at \z\ = 
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500 pc in our simulations, we still find a mass loading factor 
~ 100 at E gas = 10 Mq pc -2 , and a slope ~ —2. 

There are several differences between our simulations 
and those of Creasey et ah (2013) that could explain these 
discrepancies. Most importantly, Creasey et al. (2013) did 
not include radiative cooling below 10 4 K, whereas gas in our 
simulations can cool to 10 K. The presence of a cold phase 
in the ISM is likely to affect how outflowing gas escapes the 
disc. Creasey et al. (2013) also did not include rotation of 
the disc, or self-gravity of the gas. However, they did use a 
slightly higher resolution (1.6 pc) than we have in our simu¬ 
lations (we use a gravitational softening length of 3.1 pc for 
gas particles). Creasey et al. (2013) also injected energy from 
individual super novae, whereas our stellar feedback model 
requires that we inject energy from several supernovae si¬ 
multaneously in a single feedback event, to prevent artificial 
radiative losses (see section 2.3). 

3.3.1 Chemistry of outflowing gas 

We have seen that the mass outflow rates and velocities are 
generally not strongly affected by non-equilibrium cooling. 
However, we might expect the abundances in outflowing gas 
to be out of equilibrium, since this gas is highly dynamic. 
This would be important for comparing with observations of 
particular chemical species in outflows. For example, molec¬ 
ular outflows have been observed in extragalactic systems, 
including starbursting galaxies and Active Galactic Nuclei 
(AGN). These have been observed in emission and absorp¬ 
tion from a number of molecules, including CO, H 2 , OH and 
HCO + (e.g. Baan et al. 1989; Walter et al. 2002; Leon et al. 
2007; Sakamoto et al. 2009; Sturm et al. 2011; Emonts et al. 
2014; Geach et al. 2014). 

To investigate whether the molecular abundances are 
out of equilibrium in outflowing gas, we measured the total 
mass of H 2 in particles with a vertical velocity (perpendic¬ 
ular to the plane of the disc) greater than some velocity 
v z , i.e. Mh 2 (> v z ). We include only particles moving away 
from the mid-plane of the disc (i.e. that are outflowing), 
with z • v > 0. 

In Fig. 7 we plot Mh 2 (> v z ) as a function of v z , aver¬ 
aged over ten snapshot outputs for each simulation at inter¬ 
vals of 100 Myr. We show different metallicities and different 
radiation fields in the top and bottom panels respectively. 
The solid curves in each panel show the H 2 mass computed 
from the simulations evolved with the full chemical model 
using non-equilibrium H 2 abundances, while the dotted 
curves are calculated from the same non-equilibrium sim¬ 
ulations, but using abundances that are set to equilibrium 
in post-processing. The latter we label as ‘NonEq_Eqm’, 
to distinguish them from the ‘Eqm’ simulations that were 
evolved with cooling rates in chemical equilibrium (not 
shown). By comparing the ‘NonEq’ and ‘NonEq_Eqm’ mod¬ 
els in Fig. 7, we consider the same distribution of gas in a 
given snapshot, and the only differences are due to the H 2 
abundances being out of equilibrium. However, if we com¬ 
pared these with the ‘Eqm’ simulations, we could also see 
differences in the outflow rate of molecular hydrogen that 
are caused by the different gas distribution in that simula¬ 
tion. For example, if there had recently been a strong burst 
of star formation in one snapshot that did not occur in the 
other simulation, we would see an enhanced molecular out- 



v z (km s x ) 


Figure 7. Mass of molecular hydrogen in particles moving away 
from the mid-plane of the disc with a vertical velocity > v z , plot¬ 
ted as a function of v Zl for different metallicities (top panel) and 
for different radiation fields (bottom panel; see Table 1 ). We cal¬ 
culate the H 2 masses from simulations run using the full chemical 
model using non-equilibrium abundances (solid curves) and from 
the same simulations but with abundances set to chemical equilib¬ 
rium in post-processing (dotted curves). We average these curves 
over ten snapshots for each simulation, at intervals of 100 Myr. 
The mass of H 2 in outflowing gas (> 50kms _1 ) is generally much 
higher in non-equilibrium, e.g. by a factor of ~ 20 in the ISRF 
run. 


flow that was not caused by non-equilibrium effects. There¬ 
fore, we do not include the ‘Eqm’ simulations in Fig. 7. 

At low v z (below 5kms _1 ) the non-equilibrium H 2 
masses are generally slightly lower than in equilibrium. The 
H 2 mass at low v z is dominated by molecular clouds in the 
disc and, as we will see in section 3.4, H 2 is underabundant in 
gas that is starting to become molecular. However, at larger 
v z , where we include only gas that is outflowing, we find 
much larger H 2 masses when we use non-equilbrium abun¬ 
dances than when we assume chemical equilibrium. This gas 
was previously in molecular clouds, but has not yet had 
enough time for the H 2 to be destroyed and reach the new 
chemical equilibrium state since being ejected. For exam¬ 
ple, in the simulation at 0.1 Z 0 in the presence of the Black 
(1987) ISRF (black curves), we find 600 M 0 of molecular hy¬ 
drogen outflowing at > 50kms _1 , compared to only 30 M 0 
if we assume chemical equilibrium. Thus, non-equilibrium 
chemistry enhances the mass of outflowing H 2 by a factor 
~ 20 in this example. We see such differences in all runs 
with the Black (1987) ISRF and ten per cent of the Black 
(1987) ISRF, but the differences due to non-equilibrium 
abundances are much smaller in the presence of the Haardt 
& Madau (2001) UVB (blue curves in the bottom panel). 
Also, in the simulation without self-shielding (green curves 
in the bottom panel), the non-equilibrium and equilibrium 
H 2 masses are almost identical, and are much lower (by more 
than two orders of magnitude) than in the corresponding run 
with self-shielding, because gas does not become shielded 
from the dissociating radiation in this run. 

Fig. 7 therefore demonstrates that non-equilibrium 


MNRAS 000, 1-26 (2016) 












14 A. J. Richings and J. Schaye 



-6 -4 -2 0 2 4 6 -4 -2 0 2 4 6 -4 -2 0 2 4 6 

l°gio n H ( cm ' 3 ) 


10“ 3 


1(T 4 

10" 5 

10“ 6 

10“ 7 


PI 

O 

• i—i 

o 


a 


Figure 8. The distribution of gas in the temperature-density plane from simulations at different metallicities (top row) and for different 
radiation fields (bottom row; see Table 1 ). These simulations were run using the full non-equilibrium chemical model. In each panel we 
combine snapshot outputs taken at intervals of 100 Myr. The colour scale indicates the mass fraction of gas in each pixel. The dotted 
horizontal and vertical lines indicate the temperature and density thresholds, respectively, of our star formation prescription. Thus gas 
in the bottom right region of each panel is allowed to form stars. We find more cold, star forming gas in the simulations run at higher 
metallicity and in the presence of a weaker UV radiation field. 


chemistry can be very important for modelling molecular 
outflows, as the molecules in gas that is ejected from the 
galaxy are not instantly destroyed, but instead take time to 
evolve to a new chemical equilibrium. 

3.4 Phase structure of the ISM 

Fig. 8 shows two-dimensional histograms of the gas tem¬ 
perature and density for simulations with different metal¬ 
licities (top row) and radiation fields (bottom row), evolved 
with the full non-equilibrium chemical model. In each panel 
we stack snapshot outputs taken at intervals of 100 Myr 
to show the time-averaged distribution of gas. The colour 
scale indicates the mass fraction of gas in each pixel, so we 
see the mass-weighted distribution, rather than the volume- 
weighted distribution. 

Comparing different metallicities (top row), we see that 
increasing the metallicity allows more gas to cool to lower 
temperatures and higher densities. Reducing the strength 
of the radiation field (top centre, bottom left and bottom 
centre panels) also increases the amount of cold gas in the 
simulation, due to the reduced photoionisation heating and 
photoelectric dust heating. Note that, in the star formation 
prescription that we use, gas particles are allowed to form 


stars if they have a density nu > 1.0 cm -3 and a temper¬ 
ature T < 10 3 K, i.e. if they lie in the bottom right region 
delineated by the vertical and horizontal dotted lines in each 
panel. We thus see that there is more star forming gas in the 
galaxies at higher metallicity or in the presence of a weaker 
UV radiation field. This explains the higher star formation 
rates that we saw in section 3.2, and hence the stronger 
supernova-driven galactic winds that we saw in section 3.3, 
at higher metallicity and lower radiation field strength. 

At densities 10 -2 cm -3 <nn < 10° cm -3 , the gas at low 
metallicity follows two distinct tracks in the temperature- 
density plane, which become less distinct at higher metal¬ 
licity. We find that the high-temperature track consists of 
strongly ionised gas, which experiences strong photoheat¬ 
ing, while the low-temperature track consists of neutral gas 
that has become shielded from hydrogen-ionising radiation 
(see, for example, the first and third rows of Fig. 10, which 
show temperature-density diagrams with a colour scale indi¬ 
cating the electron abundance). At high metallicities, metal 
cooling reduces the thermal equilibrium temperature of the 
ionised gas and thus brings the two tracks closer together. 
We also see that the distinction between these two tracks is 
less clear in the presence of a weaker UV radiation field, and 
that they coincide if we neglect self-shielding. 
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Figure 9. One-dimensional mass-weighted density pdf of cold 
gas (T < 10 3 K) from simulations evolved using the full non- 
equilibrium chemical model (solid curves) and using cooling rates 
in chemical equilibrium (dashed curves) at different metallicities 
(top panel) and for different UV radiation fields (bottom panel; 
see Table 1 ). The vertical dotted lines show the minimum density 
of the cold neutral medium predicted by the model of Wolfire et 
al. (2003), as a function of metallicity and radiation field. 


In the bottom right panel of Fig. 8, we see that there is 
much less scatter in the gas temperature and density when 
we do not include self-shielding of gas from the radiation 
field, with most of the gas at densities nu > 10“ 1 cm -3 
following a very narrow region in the temperature-density 
plane. When we include self-shielding, gas particles at a par¬ 
ticular density and temperature can exhibit a wide range 
of different shielding column densities, which explains the 
larger scatter that we see in the temperature-density plane 
when self-shielding is included. 

In the simulations evolved with equilibrium cooling 
rates, we generally find similar distributions to those in 
Fig. 8 for the simulations evolved with the full non¬ 
equilibrium chemical model. However, there are some re¬ 
gions in this plane that appear different. To quantify these 
differences, we show the one-dimensional probability den¬ 
sity functions (pdfs) of gas density in Fig. 9 for simulations 
with different metallicities (top panel) and different radia¬ 
tion fields (bottom panel). Solid and dashed curves corre¬ 
spond to the full non-equilibrium chemical model and cool¬ 
ing rates in chemical equilibrium respectively. We show here 
the density distributions for cold gas, with T < 10 3 K, which 
corresponds to the temperature threshold below which gas 
particles can form stars in our star formation prescription. 
This temperature is indicated by the dotted horizontal lines 
in Fig. 8. 


At metallicities Z ^ 0.1 Z®, we see that the density dis¬ 
tribution in the top panel of Fig. 9 extends to lower densities, 
by ~ 0.5 dex, when the galaxy is evolved using the full non¬ 
equilibrium chemical model. We also see this effect when we 
consider weaker UV radiation fields at 0.1 Z®, shown by the 
red and blue curves in the bottom panel. 

For comparison, the vertical dotted lines show the min¬ 
imum density of the cold neutral medium (CNM) predicted 
by the model of Wolfire et al. (2003), if the CNM is in pres¬ 
sure equilibrium with the warm neutral medium (WNM). 
To determine this minimum density, Wolfire et al. (2003) 
calculate the thermal equilibrium temperature of the gas as 
a function of density, assuming that photoelectric heating 
from dust grains is balanced by radiative cooling from Cn 
and Oi. Using this temperature-density relation, they deter¬ 
mine the minimum pressure at which two stable ISM phases 
can exist in pressure equilibrium. This minimum pressure 
corresponds to the minimum density of the CNM, which is 
given by their equation 35: 


-31 


G'oiZ'JZ^) 


l + 3.1(G' 0 Z' d /CX 365 


(3.4) 


where Go is the strength of the UV radiation field in units 
of the Habing (1968) field, Zd is the abundance of dust and 
polyaromatic hydrocarbons (PAHs), Z g is the gas phase 
metallicity, and £t is the ionisation rate from cosmic rays 
and EUV/X-ray radiation. Primes indicate that these val¬ 
ues have been normalised to their values in the local solar 
neighbourhood, for which Wolfire et al. (2003) take Go = 1.7 
and Ct — 10 -16 s _1 . 

At solar metallicity, the minimum CNM density pre¬ 
dicted by Wolfire et al. (2003) coincides with the peak of 
the density distribution from our simulations, below which 
the distribution declines rapidly. Since this low-density tail 
is more extended in the simulations run with the full non¬ 
equilibrium chemical model, we find more cold gas below 
the minimum density of the Wolfire et al. (2003) model in 
this case than when we evolve the galaxy with equilibrium 
cooling rates. For example, at Z = Z®, 27.9 per cent of the 
cold gas mass has a density n < n m i n when we use non¬ 
equilibrium cooling rates, compared to 17.9 per cent when 
we use cooling rates in chemical equilibrium. 

Note that, in the Wolfire et al. (2003) model, the CNM 
has a maximum temperature of 243 K, whereas we consider 
cold gas with T < 10 3 K, since this corresponds to the tem¬ 
perature threshold that we use in our star formation pre¬ 
scription. If we consider the density distribution of gas with 
T < 243 K in our simulations, we find that 2.4 per cent of the 
gas mass has a density n < n m i n in the non-equilibrium run 
at solar metallicity, compared to 0.6 per cent in the equilib¬ 
rium run. Thus our simulations are not in conflict with the 
predictions of Wolfire et al. (2003). We continue to use n m i n 
as a convenient metallicity- and radiation field-dependent 
reference point in our comparisons below. 

At lower metallicities, the peak of the density distribu¬ 
tion of cold gas moves to higher densities. The minimum 
density predicted by Wolfire et al. (2003) also increases as 
the metallicity decreases, although it does not increase as 
quickly as in our simulations. We thus find less cold gas with 
n < rimin at lower metallicity. For example, at Z — 0.1 Z®, 
14.2 per cent of the cold (T < 10 3 K) gas mass has n < n m i n 
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when we use non-equilibrium cooling, compared to 6.5 per 
cent when we use cooling rates in chemical equilibrium. 

In the presence of a weaker UV radiation field (i.e. lower 
G o), the density distribution of cold gas moves to lower den¬ 
sities in our simulations, as does the minimum CNM density 
predicted by Wolfire et al. (2003). We also see that the dif¬ 
ferences between non-equilibrium and equilibrium cooling 
become more pronounced in the presence of weaker UV ra¬ 
diation fields. For example, in the simulations run in the 
presence of the Haardt V Madau (2001) UVB (blue curves 
in the bottom panel of Fig. 9), 13.2 per cent of the cold gas 
mass in the simulation run with non-equilibrium cooling has 
a density below n m in predicted by equation 3.4, compared to 
0.8 per cent in the corresponding simulation evolved using 
cooling rates in chemical equilibrium. 

To understand why non-equilibrium cooling causes 
these differences in the temperature-density distribution, it 
is useful to explore in which regions of the temperature- 
density plane the chemical abundances are out of equilib¬ 
rium. In Fig. 10 we show temperature-density diagrams with 
a colour scale indicating the mass-weighted mean electron 
abundance (first and third rows) and the mass-weighted 
mean ratio between the non-equilibrium electron abundance 
and the abundance of electrons in chemical equilibrium (sec¬ 
ond and fourth rows). Red and blue in the second and fourth 
rows indicate that the electron abundance is enhanced and 
reduced respectively, with respect to equilibrium. The top 
two rows and bottom two rows of Fig. 10 show variations in 
metallicity and radiation field respectively. 

At metallicities Z ^ 0.1 Z® the electron abundance is 
close to equilibrium throughout most of the simulation. In 
the centre panel of the second row of Fig. 10 the electron 
abundance at densities tih ~ 10 2 cm -3 is enhanced by a fac¬ 
tor of a few compared to equilibrium, and cold gas (with 
T < 10 3 K) at densities tih ~ 1 cm -3 shows electron abun¬ 
dances that are reduced by ~ 20 — 40 per cent below equi¬ 
librium. The narrow region between the warm ionised and 
warm neutral phases at T ~ 10 4 K also shows electron abun¬ 
dances that are out of equilibrium, by up to an order of 
magnitude. 

Similar trends are also seen in the presence of a weaker 
radiation field, in the left and centre panels of the bottom 
row of Fig. 10. When we neglect self-shielding, in the right 
panels of the bottom two rows of Fig. 10, the electron abun¬ 
dance is everywhere consistent with chemical equilibrium. 

At solar metallicity, in the right panels of the top two 
rows of Fig. 10, there is also a region at tih ~ 0.1 cm -3 that 
extends from T ~ 250 K to T ~ 4000 K where the electron 
abundance is enhanced by about an order of magnitude. We 
find that there is no gas in this region of the temperature- 
density plane in the simulation evolved with cooling rates 
in chemical equilibrium. In the non-equilibrium simulation, 
the enhanced electron abundance increases the collisional 
excitation rate of ions (such as Cn) by electrons and thus 
increases the radiative cooling rate, allowing it to cool be¬ 
low the thermal equilibrium temperature corresponding to 
chemical equilibrium. 

In Fig. 11 we similarly show temperature-density di¬ 
agrams with a colour scale indicating the non-equilibrium 
H 2 abundance (first and third rows) and the ratio between 
the non-equilibrium and equilibrium H 2 abundances (second 


and fourth rows), for different metallicities (top two rows) 
and different radiation fields (bottom two rows). 

At metallicities Z ^ 0.1 Z®, in the centre and right pan¬ 
els of the second row of Fig. 11, and at lower radiation field 
strengths, in the left and centre panels of the bottom row of 
Fig. 11, we see that gas with T < 10 3 K and tih ~ 1cm -3 
has a very strongly enhanced H 2 abundance, by up to six 
orders of magnitude. We find that these gas particles were 
previously in molecular clouds, with densities ~ 10 — 100 
times their current value in the previous ~ 5Myr. These 
molecular clouds then became disrupted or destroyed, and 
the gas moved to lower densities. However, since it takes a 
finite time for the molecular hydrogen to be destroyed and 
reach a new equilibrium, they retain a high molecular frac¬ 
tion, resulting in a strong overabundance of H 2 with respect 
to equilibrium. 

Similar non-equilibrium effects in the H 2 fraction were 
also found by Dobbs et al. (2008) in their simulations of 
spiral galaxies. However, it is possible that this enhance¬ 
ment in the H 2 abundance is sensitive to the resolution of 
our simulations. For example, this low-density gas may in 
reality be clumpy on scales smaller than we resolve, which 
would make it less well shielded from the photodissociating 
radiation than in our simulations. 

This region of enhanced H 2 in the temperature-density 
plane corresponds to the extended low-density tail that we 
saw in the density distributions of cold (T < 10 3 K) gas in 
Fig. 9. The enhanced H 2 abundance increases the cooling 
rate from H 2 , which allows this low-density gas to cool to 
lower temperatures. The enhanced H 2 abundances in out¬ 
flowing gas that we saw in Fig. 7 also correspond primarily 
to this region of the temperature-density plane. 

At higher densities, we also see blue regions in the sec¬ 
ond row of Fig. 11 at all metallicities, and in the left and 
centre panels of the bottom row, where the H 2 abundance 
is reduced below its equilibrium value by up to an order 
of magnitude. We find that these gas particles were previ¬ 
ously at lower densities and higher temperatures, and are 
now starting to become molecular. 

Pelupessy V Papadopoulos (2009) also explored non¬ 
equilibrium H 2 chemistry in their hydrodynamic simula¬ 
tions of isolated galaxies, with metallicities up to Z®. They 
found H 2 fractions that were often far out of equilibrium in 
their simulations, in agreement with our results. In contrast, 
Krumholz V Gnedin (2011) compared the time-dependent 
H 2 model from Gnedin V Kravtsov (2011) to the equilibrium 
H 2 model of Krumholz et al. (2008, 2009) and McKee V 
Krumholz (2010), implemented in cosmological simulations, 
and they found good agreement between these two models at 
metallicities > 0.01 Z®. They therefore concluded that their 
equilibrium treatment of H 2 is sufficient above one per cent 
solar metallicity. However, their simulations used a lower 
resolution than we use. For example, the maximum reso¬ 
lution in their zoom-in cosmological simulations was 65 pc, 
whereas our simulations use a gravitational softening length 
of 3.1 pc. This difference in resolution may explain why we 
find non-equilibrium chemistry to be more important than 
Krumholz & Gnedin (2011). 
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Figure 10. Temperature-density diagrams from simulations run with the full non-equilibrium chemical model at different metallicities 
(top two rows) and for different radiation fields (bottom two rows; see Table 1). The colour scale indicates the mass-weighted mean 
non-equilibrium electron abundance (first and third rows) and the mass-weighted mean ratio between the non-equilibrium electron 
abundance and the electron abundance in equilibrium (second and fourth rows). We have averaged over snapshot outputs taken at 
intervals of 100 Myr. Red regions in the second and fourth rows show where electrons are enhanced with respect to equilibrium, and blue 
regions show where it is suppressed. 


4 OBSERVABLE LINE EMISSION 

To investigate the effects of metallicity, radiation field and 
non-equilibrium chemistry on observable diagnostics, we 
computed line emission maps for Cn and CO. We used 
the publicly available Monte-Carlo radiative transfer code 
radmc-3d 5 (version 0.38) to compute the line emission from 
these species by post-processing the output from our simu¬ 
lations. radmc-3d includes thermal emission from dust and 
line emission from user-specified species and transitions. We 
also include anisotropic scattering of continuum and line 
emission by dust grains, although in the current version of 
radmc-3d scattering of line emission does not include the 
corresponding doppler shift due to the relative motion of the 
dust, it only changes the direction of the radiation. 


5 http://www.ita.uni-heidelberg.de/~dullemond/software/ 
radmc-3d/ 


We include two populations of dust grains, graphite and 
silicate, using opacities from the calculations of Martin & 
Whittet (1990), who used the power-law grain size distri¬ 
bution of Mathis et al. (1977). We use a dust-to-gas mass 
ratio of 2.4 x 10 -3 Z/Z© and 4.0 x 10 -3 Z/Z© for the graphite 
and silicate grains respectively, which we take from the ‘ISM’ 
grain abundances used in the photoionisation code cloudy 6 
version 13.01 (Ferland et al. 2013). 

The line emission from a given species depends on its 
level populations. However, unless we assume that these 
level populations are in Local Thermodynamic Equilibrium 
(LTE), they will also depend on the radiation field, includ¬ 
ing the emission lines themselves. A full, self-consistent non- 
LTE treatment can therefore be computationally expensive, 
as line emission from one gas cell can influence the level pop¬ 
ulations of its neighbours. Nevertheless, non-LTE effects can 

6 http://nublado.org/ 
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Figure 11. As Fig. 10, but with the colour scale now indicating the mass-weighted mean non-equilibrium H 2 abundance (first and third 
rows) and the mass-weighted mean ratio between the non-equilibrium H 2 abundance and the H 2 abundance in equilibrium (second and 
fourth rows). We see that the H 2 abundance at T < 10 3 K and wh ~ 1cm -3 is strongly enhanced, by up to six orders of magnitude, at 
Z ^ 0.1 Zq. The resulting increase in H 2 cooling in this region explains why we see more cold gas at low densities in the one-dimensional 
density pdfs in Fig. 9 when we use non-equilibrium cooling compared to equilibrium cooling. 


be important. Duarte-Cabral et al. (2015) created synthetic 
maps of CO emission from their hydrodynamic simulations 
of spiral galaxies, and they compared maps created with and 
without assuming LTE. They found that, while the morphol¬ 
ogy of the CO emission was unaffected, the CO line intensity 
was generally overestimated when assuming LTE. 

radmc-3d includes a number of approximate methods 
to include non-LTE effects. We use the Large Velocity Gra¬ 
dient (LVG) method, which is also known as the Sobolev 
approximation (Sobolev 1957). This method assumes that, 
after propagating for some distance, an emission line will be 
sufficiently Doppler shifted, due to motions of the gas, to 
propagate freely. We can thus calculate an escape probabil¬ 
ity for emitted photons based on the velocity gradient, which 
allows us to calculate the level populations from local quan¬ 
tities. A full description of the LVG method as implemented 
in radmc-3d can be found in Shetty et al. (2011). 

With the LVG method, the non-LTE line emission from 
species i depends on the density of i and the species j that 


can collisionally excite it, along with the gas temperature, 
gas velocity and the radiation field. Apart from the radia¬ 
tion field, we compute these quantities, and the densities of 
graphite and silicate grains, on a 4.0 kpc x 4.0 kpc x 4.0 kpc 
Cartesian grid with a resolution of 10 pc per cell. We calcu¬ 
late these quantities from the SPH particles in the snapshot 
outputs using SPH interpolation with the same Wendland 
C2 kernel with 100 SPH neighbours as was used in the sim¬ 
ulations. We then compute the line and thermal dust emis¬ 
sion, viewing the disc face-on, in 80 wavelength bins centred 
on the line and covering a velocity range of ±40kms -1 . Fi¬ 
nally, we repeat this without the emission line to create a 
map of the thermal dust emission only, and we subtract this 
from the full map to obtain the continuum-subtracted line 
emission. 

It is important to note that the Cn and CO emission 
presented in this section may be affected by the resolution of 
our simulations. In particular, both Cn and CO have a crit¬ 
ical density of ~ 10 3 cm -3 . However, we saw in Fig. 8 that 
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Figure 12. Cn 158 /im line emission from simulations evolved 
with the full non-equilibrium chemical model at different metal- 
licities (top row) and for different radiation fields (bottom row; see 
Table 1 ). Each map shows the central region of the galaxy, 4.0 kpc 
across, at 500 Myr, viewing the disc face on. Cn emission increases 
with increasing metallicity and increasing radiation field. 

structures with such high densities are not well resolved in 
our simulations. Therefore, we might underestimate the Cn 
and CO emission from dense, unresolved structures. Addi¬ 
tionally, high-resolution simulations of dense clouds find that 
most CO is concentrated in compact (~ lpc), high density 
10 3 cm -3 ) clumps and filaments (e.g. Glover & Clark 
2012 ). 

4.1 CII fine-structure line emission 

The Cn fine-structure line at 158 /xm is an important coolant 
in neutral, atomic gas (e.g. Richings et al. 2014a). This line 
is therefore commonly used as an observational tracer of 
the cooling properties and physical conditions of the neutral 
phases of the ISM (e.g. Malhotra et al. 2001; Brauher et 
al. 2008; Gracfa-Carpio et al. 2011; Kennicutt et al. 2011; 
Beirao et al. 2012; Croxall et al. 2012). 

We calculate the emission from the Cn 158 /im line using 
radmc-3d, as described above, with atomic data for the Cn 
ion taken from the lamda database 7 (Schoier et al. 2005). 
These include excitation rates from collisions with ortho- 
and para-H2 (Lique et al. 2013; Wiesenfeld & Goldsmith 
2014), for which we assume an ortho-to-para ratio of 3:1, 
neutral Hi (Barinovs et al. 2005), and electrons (Wilson & 
Bell 2002). 

Fig. 12 shows Cn line emission maps of the central re¬ 
gion of each galaxy, 4.0 kpc across, at 500 Myr, viewing the 
disc face on. The top row shows variations in metallicity, 
while the bottom row shows different UV radiation fields. 
These maps were computed from simulations evolved with 
the full chemical model, using the non-equilibrium abun¬ 
dances of Cn and the collision species (H 2 , Hi and electrons). 

Comparing the galaxies at different metallicities Z, we 
see that the total Cn line emission increases approximately 

7 http://home.strw.leidenuniv.nl/~moldata/ 


linearly with Z. We might have expected such a linear rela¬ 
tionship as increasing Z will increase the amount of carbon 
in the gas. However, this simple picture is complicated by 
the fact that the densities and temperatures of the neutral 
gas are also affected by the metallicity, as we saw in Fig. 8, 
which will affect the emissivities per Cn ion. A better way to 
explain this trend of Cn emission with metallicity is to note 
that, if the neutral gas is in thermal equilibrium, the cool¬ 
ing rate of the gas (which comes primarily from Cn and Oi 
line emission) will balance the heating rate (which is mainly 
from photoelectric heating from dust grains). The photoelec¬ 
tric heating rate (Bakes V Tielens 1994; Wolfire et al. 1995) 
scales with the abundance of dust grains, which we assume 
scales linearly with Z. Therefore, at higher metallicity, the 
photoelectric heating rate increases, and thus we need more 
Cn emission to achieve thermal balance. 

We also see very different morphologies of Cn emission 
at different metallicities. At Z ^ 0.1 Z©, we see Cn emission 
from dense clumps that are loosely arranged into spiral arms. 
In contrast, at Z = 0.01 Z© we see an almost undisturbed 
disc, with Cn emission peaking in the centre. This reflects 
the different morphologies of the total gas distribution that 
we saw in Fig. 1, which are due to the disruption of the disc 
by stellar feedback and the ability of the gas to cool to lower 
temperatures and higher densities at higher metallicities, as 
discussed in section 3.2. 

Comparing galaxies with different UV radiation fields 
at ten per cent solar metallicity (top centre, bottom left and 
bottom centre panels of Fig. 12), we see that the total Cn 
emission increases as the strength of the radiation field in¬ 
creases. This can also be understood as the photoelectric 
heating rate is higher in the presence of a stronger UV field, 
and thus more Cn emission is needed to balance the heat¬ 
ing rate. However, the photoelectric heating rate does not 
increase linearly with radiation field strength, because dust 
grains become positively charged in the presence of a strong 
UV field, which reduces the photoelectric heating efficiency. 

In the run without self-shielding (bottom right panel), 
we find much stronger Cn emission, by nearly an order of 
magnitude compared to the corresponding simulation that 
includes self-shielding (bottom centre panel). This is be¬ 
cause the dense gas is much warmer when self-shielding is 
ignored. For example, at nu = 100 cm -3 , T ~ 500 K when 
self-shielding is not included, compared to T ~ 50 — 100 K 
when it is included (see Fig. 8). Another way to see this 
is that, without self-shielding, the hydrogen-ionising radia¬ 
tion persists in the dense gas, so we have additional heating 
from photoionisation of hydrogen. We therefore need more 
Cn emission to achieve thermal balance. 

One caveat to note here is that, since we use a uni¬ 
form UV radiation field, and we do not model local sources 
of radiation such as young stars, our simulations do not in¬ 
clude Hii regions and PDRs. However, these are often strong 
sources of Cn emission. Therefore, it is possible that we are 
underestimating the Cn emission. 

We saw in section 3.4 that evolving a galaxy using cool¬ 
ing rates in non-equilibrium can affect the distribution of gas 
densities and temperatures, compared to using cooling rates 
in chemical equilibrium. We also found that some chemi¬ 
cal species (particularly the free electrons and H 2 ) can be 
far out of equilibrium at certain densities and temperatures. 
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Figure 13. Ratio of average line intensity assuming equilibrium 
abundances, /Eqm? to that using non-equilibrium abundances, 
/NonEqj plotted against time, for different metallicities (top row) 
and different UV radiation fields (bottom row; see Table 1). 
We calculated /E qm from simulations evolved with the full non- 
equilibrium chemical model with abundances set to equilibrium 
in post-processing (NonEq_Eqm; solid curves), and from simula¬ 
tions evolved in chemical equilibrium (Eqm; dashed curves). We 
show the Cn 158 /im line (black curves) and the CO J = 1 — 0 
line (red curves). 


These two effects could potentially have an impact on the 
observable line emission from individual species. 

To investigate the impact of non-equilibrium chemistry 
on the line emission maps, we also computed these from 
equilibrium abundances in two ways. Firstly, we used the 
simulations evolved with the full non-equilibrium chemical 
model and set the chemical abundances of each gas par¬ 
ticle to chemical equilibrium (‘NonEq_Eqm’). This shows 
how non-equilibrium abundances affect the line emission. 
Secondly, we used the simulations that were evolved using 
cooling rates in chemical equilibrium (‘Eqm’). This shows 
how the different distributions of gas density and temper¬ 
ature, due to using non-equilibrium cooling, affect the line 
emission. 

Fig. 13 shows the ratio of the average line inten¬ 
sity assuming equilibrium abundances to that using non¬ 
equilibrium abundances, plotted against time. Solid and 
dashed curves show ‘NonEq_Eqm’ and ‘Eqm’ respectively. 
Cn line emission is shown by the black curves, while the 
red curves show CO line emission, which we will discuss in 
section 4.2. The top row shows simulations run at different 
metallicities, and the bottom row shows for different UV 
radiation fields. 

The black solid curves, for Cn ‘NonEq_Eqm’, are all 
very close to unity in all panels. This tells us that the abun¬ 
dances of Cn and those species that collisionally excite it 
are close to equilibrium in Cii-emitting gas for the non¬ 
equilibrium runs. The black dashed curves, from the sim¬ 
ulations evolved in chemical equilibrium, also remain close 
to unity, except in runs at ten per cent solar and solar metal- 
licity in the presence of the Black (1987) ISRF (top centre 
and top right panels), where the equilibrium emission is ~ 50 
per cent higher. 
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Figure 14. As Fig. 12, but for CO J = 1 — 0 line emission at 
2.6 mm. CO emission increases with increasing metallicity and 
decreasing UV radiation field. 


4.2 CO line emission 

Molecular hydrogen is difficult to observe directly in cold, 
molecular gas, because the lowest rovibrational levels of the 
H 2 molecule are difficult to excite at the temperatures typi¬ 
cal of molecular clouds (~ 10 K). For example, the lowest ro¬ 
tational transition in the ground vibrational state, 0 — 0 5(0), 
has an excitation energy E/k& s= 510 K. However, one of 
the next most abundant molecules after H 2 is CO, which 
can be observed at much lower temperatures than H 2 . For 
example, the lowest rotational transition of CO, J — 1 — 0, 
has an excitation energy of 5.53 K. CO emission is there¬ 
fore commonly used to map cold molecular gas (e.g. Heifer 
et al. 2003; Kuno et al. 2007; Bolatto et al. 2008; Leroy et 
al. 2009). However, to determine the H 2 content from CO 
emission alone, we need to know the conversion factor, Vco, 
between CO emission and H 2 column density (see Bolatto et 
al. 2013 for a recent review). The Xco factor is commonly 
defined as: 


Tco = 7 ^ cm 2 (Kkms 1 ) \ (4.1) 

where Vh 2 is the H 2 column density and Ico is the velocity- 
integrated intensity of the CO J — 1 — 0 line. 

In this section we present the emission from the CO 
J — 1 — 0 line, at a wavelength of 2.6 mm, in our simu¬ 
lations, computed using radmc-3d. We use molecular CO 
data from the LAMDA database, including collisional exci¬ 
tation by ortho- and para-H 2 (Yang et al. 2010), for which 
we assume an ortho-to-para ratio of 3:1. 

We show velocity-integrated line emission maps for the 
J — 1 — 0 line in Fig. 14 for different metallicities (top 
row) and different radiation fields (bottom row). In the 
galaxy with the lowest metallicity (0.01 Z Q ; top left panel 
of Fig. 14), and in the run without self-shielding (bottom 
right panel), the CO intensity is very low, and would not be 
detectable. For comparison, for the CO observations of Local 
Group galaxies in Leroy et al. (2011), the 3 a CO intensity 
threshold for the Small Magellanic Cloud is 0.25Kkms -1 . 

We see very weak CO emission in the lowest-metallicity 
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run because there is less carbon and oxygen available to form 
CO. Additionally, there is less dust shielding at lower metal¬ 
licity, which is needed to prevent the destruction of CO by 
photodissociation. Similarly, in the simulations run without 
self-shielding, we see little CO emission because photode¬ 
struction of CO is always efficient. 

In the remaining galaxies, we see that the CO emission 
is concentrated in dense clumps along the spiral arms. Com¬ 
paring panels in the top row, the CO intensity increases with 
increasing metallicity, while in the bottom left and bottom 
centre panels, the CO intensity increases with decreasing 
radiation field strength. 

The ratio of the average CO line intensity assuming 
equilibrium to that using non-equilibrium abundances is 
shown by the red curves in Fig. 13. Non-equilibrium chem¬ 
istry has a greater effect on line emission from CO than 
for Cn. For example, in the simulations run in the pres¬ 
ence of ten per cent of the Black (1987) ISRF at ten per 
cent solar metallicity (lowISRF, bottom left panel), the av¬ 
erage CO line intensity computed in equilibrium is higher 
by a factor of ~ 4 than in non-equilibrium. We also see 
very large relative differences between equilibrium and non¬ 
equilibrium CO emission at the lowest metallicity (0.01 Z©; 
top left panel), although we saw in Fig. 14 that the CO emis¬ 
sion is very weak in this run, and would not be detectable 
in typical CO surveys. 

These non-equilibrium trends in CO emission are not 
always in the same direction. For example, at solar metallic¬ 
ity (top right panel of Fig. 13), the average CO line intensity 
calculated by setting the abundances to equilibrium in post¬ 
processing (‘NonEq_Eqm’, solid red curve) is lower than in 
non-equilibrium by a factor of ~ 2. This is opposite to the 
trend that we saw in the lowISRF run. We find that there 
are two competing non-equilibrium effects that act on the 
CO abundance, and hence on the CO intensity. Firstly, gas 
that is forming into molecular clouds takes a finite time to 
form CO. Such gas is underabundant in CO with respect to 
equilibrium. Secondly, when an existing molecular cloud is 
disrupted or destroyed, it takes a finite time for the CO in 
this gas to be destroyed. Hence, CO is overabundant in such 
gas. 

Therefore, the effect of non-equilibrium chemistry on 
CO emission is not a simple one, as there are two competing 
effects, which depend crucially on the thermal history of the 
gas. Thus, non-equilibrium abundances can either increase 
or decrease the CO intensity. 

In Fig. 15 we plot the H 2 column density of each pixel 
versus the velocity-integrated intensity of the CO J = 1 — 0 
line. Each pixel in the emission line maps is 10 pc across, and 
the maps cover the central 4kpc of the disc. We include ten 
snapshot outputs from each simulation, taken at intervals 
of 100 Myr. Many of the pixels have very low CO intensi¬ 
ties that would be undetectable in a realistic survey of CO 
emission in extragalactic sources. We therefore only include 
pixels with a CO intensity ico > 0.25 Kkms -1 , which cor¬ 
responds to the 3 a intensity threshold used by Leroy et al. 
(2011) for the Small Magellanic Cloud. 

We show the Ah 2 — ico relation from simulations run 
with the full non-equilibrium chemical model for different 
metallicities in the top row, and for different UV radiation 
fields in the bottom row. We do not include the simulations 
run at one per cent solar metallicity or run without self- 



Figure 15. Relation between H 2 column density, Ah 2 , and 
velocity-integrated intensity of the CO J = 1 — 0 line, /co? 
measured in pixels 10 pc across that span the central 4kpc 
of the galaxy, viewing the disc face-on. We show simulations 
run with the full non-equilibrium chemical model at different 
metallicities (top row) and for different radiation fields (bot¬ 
tom row; see Table 1 ). We include ten snapshots from each 
simulation, taken at intervals of 100 Myr. Black lines show 
the linear Ah2 — ico relation obtained using the mean Aco 
factor from the given simulation, averaged over pixels with 
Ico > 0.25Kkms _1 . We express this mean conversion factor 
as Aco, 20 = A C o/(10 20 cm -2 (Kkms -1 ) -1 ). 


shielding, as these do not show detectable CO emission. The 
black line in each panel shows the linear relation between 
iVn 2 and Ico using the mean Aco factor measured in the 
simulation, averaged over pixels with Ico > 0.25Kkms -1 . 

In Fig. 15, we see that the mean Aco factor in¬ 
creases with increasing radiation field strength, Go. For 
comparison, measurements of the Aco factor in molecu¬ 
lar clouds in the Milky Way, using various different meth¬ 
ods (Virial mass, CO isotopologues, dust extinction or 
emission, diffuse gamma-ray emission), find Aco ,20 = 

Aco/(10 20 cm -2 (Kkms -1 ) -1 ) ~ 2 — 4 (e.g. Bolatto et al. 
2013 and references therein). This trend with Go may arise 
because, in our simulations with different radiation fields 
(at 0.1 Z 0 ), the dust extinction only extends up to A v ~ 1 . 
CO is therefore not fully shielded from photodissociation by 
dust in these simulations, whereas the molecular hydrogen 
can still become fully shielded, by self-shielding. Therefore, 
as Go increases, the CO intensity is suppressed more than 
the H 2 column density, and thus Aco increases. 

In Fig. 15, there is almost no dependence of Aco on 
metallicity. This lack of dependence on metallicity may seem 
surprising, as observations find that Aco decreases with in¬ 
creasing metallicity (e.g. Israel 1997; Leroy et al. 2011). How¬ 
ever, this apparent discrepancy is likely because our simu¬ 
lations at ten per cent solar metallicity only probe regions 
with dust extinctions A v < 1 , while the solar metallicity 
runs extend up to A v ~ 10. 

To understand the trends of the Aco factor at different 
A v , we can look at the model of Feldmann et al. (2012). 
They construct a model for Aco as a function of metal- 
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Figure 16. Best-fit power law relation between H 2 column den¬ 
sity, 7 Vh 2 5 and velocity-integrated intensity of the CO J = 1 — 0 
line, JcO) at different metallicities (top row) and for different ra¬ 
diation fields (bottom row; see Table 1 ). We use maps of Xh 2 
and Iqq computed from simulations evolved with the full non- 
equilibrium chemical model using non-equilibrium abundances 
(NonEq; solid curves) and with abundances set to equilibrium in 
post-processing (NonEq_Eqm; dashed curves), and from simula¬ 
tions evolved in chemical equilibrium (Eqm; dot-dashed curves). 
Dotted curves indicate lines of constant JcO) as indicated in 
the top right panel. Each relation was fit to pixels with Ico > 
0.25 Kkms- 1 . 

licity, Z, radiation field, U, and dust extinction, A v . They 
use H 2 and CO abundances from small-scale MHD simula¬ 
tions of the turbulent ISM from Glover X Mac Low (2011), 
which included a treatment for the non-equilibrium chem¬ 
istry of H 2 and CO, although Feldmann et al. ( 2012 ) assume 
photodissociation equilibrium to determine the dependence 
of the CO abundance on U, as most of the simulations of 
Glover X Mac Low (2011) were run for only one radiation 
field. Feldmann et al. ( 2012 ) compute the CO line emission 
based on an escape probability formalism, assuming that the 
level populations of the CO molecule are in LTE and using 
an assumption for the CO line width (either a constant line 
width or a virial scaling). 

Feldmann et al. (2012) show the dependence of Xco on 
A v in their model for various Z and U (see their fig. 2 ). They 
show that Xco reaches a minimum at A v ~ 2 — 10, which 
approximately corresponds to where the CO line becomes 
optically thick. At A v below this minimum (where the CO 
line is optically thin), they show that, at fixed A v , Xco 
increases with U (in agreement with our simulations), but 
is independent of Z. 

In our simulations at solar metallicity, we probe regions 
up to A v ~ 10. However, all our simulations at lower metal¬ 
licities only extend up to A v ~ 1 . Therefore, it is likely that 
we do not recover the observed trends of Xco with metal¬ 
licity because we do not cover a range of metallicities in the 
optically thick regime. 

The Xco conversion factor between CO intensity and 
H 2 column density, as defined in equation 4.1, may also be 
affected by non-equilibrium chemistry, which affects both 
the CO emission (as seen in Fig. 13) and the abundance of 


Table 3. Mean Xco factor, averaged over pixels with a velocity- 
integrated CO line intensity Ico > 0.25Kkms -1 . 


Simulation 

Xco/10 20 cm 2 (Kkms 
NonEq NonEq Eqm 

i)-i 

Eqm 

ref 

10.2 

23.9 

23.0 

hiZ 

9.0 

23.4 

13.6 

lowISRF 

4.7 

5.9 

6.6 

UVB 

3.2 

3.5 

3.3 


H 2 (as seen in Fig. 11). Fig. 16 compares the best-fit power- 
law relation between 7 Vh 2 and Ico from the non-equilibrium 
simulations (‘NonEq’; solid curves) to those from the same 
simulations with abundances set to equilibrium in post¬ 
processing (‘NonEq_Eqm’; dashed curves), and from sim¬ 
ulations evolved in chemical equilibrium (‘Eqm’; dot-dashed 
curves). Dotted curves indicate lines of constant Xco- The 
top and bottom rows show different metallicities and UV ra¬ 
diation fields respectively. Each curve was fit to pixels with 
Ico > 0.25Kkms -1 . 

The best-fit relations in Fig. 16 have power-law slopes 
of rsj 0.4 — 0.6. This is flatter than the linear relation that we 
would expect for a constant Xco factor, and indicates that 
Xco decreases as Ico increases. 

The 7 Vh 2 — Ico relations calculated using non¬ 
equilibrium abundances (solid curves) are lower than those 
calculated in equilibrium. Table 3 summarises the mean 
Xco factors from each simulation, averaged over pixels with 
Ico > 0.25Kkms -1 . 

In the presence of the Black (1987) ISRF (ref and hiZ), 
we find that the mean Xco factor decreases by a factor 
~ 2.3 when we use non-equilibrium abundances. At lower 
radiation field strengths (lowISRF and UVB), the effect of 
non-equilibrium chemistry on the mean Xco factor is much 
smaller (e.g. 10 per cent in the UVB runs). Nevertheless, 
we saw in Fig. 13 that the mean CO intensity for weaker 
radiation fields is lower, by a factor of ~ 4, when we use 
non-equilibrium abundances than in equilibrium. We also 
find fewer pixels with Ico above the detection threshold of 
0.25Kkms -1 when we use non-equilibrium abundances in 
these examples. 

If we find fewer pixels with detectable CO emission 
when we use non-equilibrium abundances, we might also 
expect this to have an impact on the fraction of CO-dark 
molecular gas, i.e. molecular hydrogen that is not traced by 
CO emission (e.g. Tielens X Hollenbach 1985; van Dishoeck 
X Black 1988; Wolfire et al. 2010; Smith et al. 2014). We 
calculated the CO dark fraction in the simulations ref, hiZ, 
lowISRF and UVB. This is the mass of H 2 in pixels with 
Ico < 0.25Kkms -1 divided by the mass of H 2 in all pixels 
of the CO map. For each simulation, we averaged this CO 
dark fraction over ten snapshots taken at intervals of 100 
Myr. When we use H 2 abundances and CO maps in non¬ 
equilibrium, the CO dark fraction is 0.76 — 0.86 in these 
four simulations, with no strong dependence on metallicity 
or UV radiation field. However, when we set the abundances 
to chemical equilibrium, the CO dark fraction in the UVB 
and lowISRF runs decreases to 0.63 and 0.59, respectively. 
In the ref and hiZ runs, the effect of non-equilibrium chem- 
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istry is weaker, as the CO dark fraction only changes by ten 
per cent when we use equilibrium abundances. 

Wo hire et al. (2010) used theoretical models of molec¬ 
ular clouds to determine the CO dark fraction for different 
cloud masses and UV radiation fields. They found an almost 
constant CO dark fraction of 0.3, although they only con¬ 
sidered individual clouds, whereas we compute the CO dark 
fraction over the disc of the galaxy, including diffuse gas. 
Smith et al. (2014) simulated four Milky Way-type galaxies 
with different gas surface densities and UV radiation fields. 
They found CO dark fractions of 0.4 — 0.85. The highest 
values that they found are consistent with our simulations, 
although they also found that increasing the strength of the 
UV radiation field by a factor of 10 increases the CO dark 
fraction by a factor of « 1.5 — 2. We do not find any strong 
dependence on UV radiation field in our simulations. 


5 CONCLUSIONS 

We have run a series of hydrodynamic simulations of isolated 
galaxies with a virial mass M 2 oo,crit = 10 11 M© and stellar 
mass 10 9 M©. The models use a resolution of 750 M© per 
gas particle and a gravitational force resolution of 3.1 pc, and 
were run with a modified version of the SPH code gadget3. 
We included a treatment for the full non-equilibrium chem¬ 
ical evolution of ions and molecules (157 species in total), 
along with gas cooling rates computed self-consistently from 
these non-equilibrium abundances (Richings et al. 2014a, b), 
and we compared these to simulations evolved using cooling 
rates in chemical equilibrium. 

Our simulations were run at a fixed metallicity and in 
the presence of a uniform UV radiation field, with a local 
prescription for self-shielding by gas and dust. We covered 
a wide range of metallicities (0.01 Z©, 0.1 Z© and Z©), and 
different UV radiation fields that span nearly three orders 
of magnitude in Hi photoionisation rate (the Black (1987) 
ISRF, ten per cent of the Black (1987) ISRF, and the redshift 
zero UVB of Haardt & Madau (2001); see Table 1), and we 
also repeated the runs with the Haardt & Madau (2001) 
UVB without self-shielding. 

Our goal was to investigate the effects of metallicity, ra¬ 
diation field and non-equilibrium chemistry, which have all 
been demonstrated to affect gas cooling rates, in simulations 
of galaxy evolution. There are two aspects to the impact of 
these effects. Firstly, how the changes in gas cooling rates 
affect the evolution of the galaxy, and secondly, how observ¬ 
able tracers of individual chemical species are affected. 

Our main results are as follows: 

(i) In simulations at higher metallicity, and for weaker 
UV radiation fields, gas can more easily cool to a cold 
(T rsj 100 K), star forming ISM phase, due to increased 
metal-line cooling and reduced UV heating respectively. We 
thus find higher star formation rates in these cases, by two 
orders of magnitude and a factor ~ 3 for the different metal¬ 
licities and different radiation fields, respectively, that we 
consider here (Fig. 4). In particular, the gas surface density 
threshold below which star formation is cut off decreases 
with increasing metallicity and decreasing UV radiation field 
(Fig. 3), as predicted by Schaye (2004). 

(ii) We find higher mass outflow rates at higher metallic¬ 
ity (by two orders of magnitude) and for weaker radiation 


fields (by a factor ~ 3), due to the higher star formation 
rates (Fig. 5). However, the average mass loading factor (i.e. 
the ratio of outflow to star formation rates), measured at 
0 . 2 R 2 oo,crit = 19kpc, is ~ 10, regardless of metallicity or 
radiation field. 

(iii) The mass loading factor, /3, measured 1 kpc above 
and below the disc decreases with increasing gas surface den¬ 
sity, Egas, following approximately a power-law, ft oc E - a g 
(Fig. 6). At fixed E gas , the mass loading factor increases 
with decreasing metallicity and increasing radiation field 
strength. 

(iv) Non-equilibrium cooling does not strongly affect the 
total star formation rate of the galaxy (Fig. 4). The ini¬ 
tial rise in mass outflow rate in the first ~ 200 Myr of the 
simulation is sometimes more gradual when we use equilib¬ 
rium cooling, compared to the non-equilibrium runs (Fig. 5). 
However, apart from this initial difference, non-equilibrium 
cooling does not strongly affect the outflow properties. 

(v) Non-equilibrium chemistry does have a large effect on 
the chemical make-up of outflowing gas (Fig. 7). For exam¬ 
ple, in our reference run (black curves in Fig. 7), we find 
on average 600 M© of H 2 outflowing with a vertical veloc¬ 
ity of > 50kms -1 if we use non-equilibrium abundances, 
compared to 30 M© if we assume chemical equilibrium. Non¬ 
equilibrium chemistry therefore enhances the mass of out¬ 
flowing H 2 by a factor ~ 20 in this example. This has im¬ 
portant implications for modelling molecular outflows in hy¬ 
drodynamic simulations of galaxies. 

(vi) We investigated where in the temperature-density 
plane molecular hydrogen is out of equilibrium (Fig. 11). H 2 
can be enhanced, by up to six orders of magnitude, in gas 
that was previously in molecular clouds at higher densities 
but has since been disrupted. We also find regions, around 
n h ~ 10 — 100cm -3 , where H 2 is underabundant, by up to 
an order of magnitude. This is due to gas that is starting to 
form molecular clouds, but has not yet had enough time to 
fully form H 2 . 

(vii) Using the publicly available Monte-Carlo radiative 
transfer code radmc-3d 8 , we performed radiative transfer 
calculations on our simulations in post-processing to com¬ 
pute the line emission from Cn and CO. Cn emission from 
the 158 /im line is stronger at higher metallicity and for 
stronger radiation fields (Fig. 12), while CO emission from 
the J = 1 — 0 line is stronger at higher metallicity and for 
weaker radiation fields (Fig. 14). 

(viii) Cn emission is generally unaffected by non¬ 
equilibrium chemistry, whereas CO emission is affected by a 
factor of ~ 2 — 4 (Fig. 13). However, the CO emission can 
be either higher or lower in non-equilibrium, since, similarly 
to H 2 , the CO abundance can be either enhanced or sup¬ 
pressed. This also affects the mean Xco conversion factor 
between CO line intensity and H 2 column density (equa¬ 
tion 4.1) that we measure in the simulations, by up to a 
factor ~ 2.3 (Table. 3). 

(ix) Non-equilibrium chemistry also affects the fraction 
of CO-dark molecular gas (e.g. Tielens & Hollenbach 1985; 
van Dishoeck & Black 1988; Wolfire et al. 2010; Smith et 
al. 2014), i.e. the fraction of molecular hydrogen that is not 


8 http://www.ita.uni-heidelberg.de/~dullemond/software/ 
radmc-3d/ 
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traced by observable CO emission. For example, in our c low- 
ISRF’ run, we find 86 per cent of the H 2 mass in the central 
4 kpc of the disc lies in pixels with Ico < 0.25 K km s -1 if we 
use non-equilibrium abundances, compared to 59 per cent if 
we assume chemical equilibrium. 

To summarise, we have demonstrated that metallic- 
ity and UV radiation affect the global properties of galax¬ 
ies, such as their star formation rates, and the observable 
signatures of individual chemical species. In contrast, non¬ 
equilibrium chemistry and cooling generally do not strongly 
affect the global properties of galaxies, although they do 
affect the observable diagnostics, particularly in molecular 
gas. 

To investigate how sensitive these results are to the res¬ 
olution of the simulations, we also repeated the ref simu¬ 
lation twice, with a mass resolution a factor of four higher 
and lower than our fiducial resolution of 750 M©. We find 
that our results are mostly unaffected by resolution over 
the range that we test here. In particular, the star forma¬ 
tion rates (Fig. 4) are similar at the three different resolu¬ 
tions, although the low resolution run shows a smoother star 
formation history with no strong bursts. The mass loading 
factors of the outflows (Fig. 6) are also similar at differ¬ 
ent resolutions. The mass of outflowing molecular hydrogen 
(Fig. 7) is similar at our fiducial and high resolutions, al¬ 
though the low resolution run shows less difference between 
non-equilibrium and equilibrium in the molecular outflows. 

However, there are several important caveats that we 
need to highlight that may mean that we have underesti¬ 
mated the importance of non-equilibrium chemistry in these 
simulations. Firstly, we simulate isolated galaxies that do 
not include cosmological processes such as galaxy merg¬ 
ers and accretion of gas onto the galaxy. We might expect 
that such processes could enhance the non-equilibrium ef¬ 
fects. For example, during and immediately after a merger, 
the ISM will evolve rapidly as it re-distributes its gas in 
the temperature-density plane. This could potentially drive 
chemical abundances further out of equilibrium. To investi¬ 
gate the importance of such processes for the chemistry, we 
could repeat this study using cosmological zoomed simula¬ 
tions of individual galaxy haloes. 

Secondly, turbulence in the ISM can also drive chemi¬ 
cal abundances out of equilibrium. For example, Gray et al. 
(2015) recently presented a series of high-resolution simula¬ 
tions of the turbulent ISM. They showed that the steady- 
state ion abundances at the end of their simulations can be 
out of equilibrium by several orders of magnitude at high 
Mach numbers. If, as expected, we do not fully resolve such 
small-scale turbulence in our simulations, then we are likely 
to underestimate the importance of non-equilibrium chem¬ 
istry. 

Thirdly, chemical abundances can also be driven out 
of equilibrium by a fluctuating radiation field. For exam¬ 
ple, Oppenheimer & Schaye (2013b) demonstrated that the 
abundances of metal ions in the circumgalactic medium can 
be affected by the presence of an AGN even after the AGN 
has turned off, as it takes a long time for the ions to recom¬ 
bine. In our simulations, we apply a local prescription for 
self-shielding of the radiation field by gas and dust, which 
does vary with position and time. However, we apply the 
self-shielding to a constant, uniform radiation field. In real¬ 


ity, the UV radiation from young stars will fluctuate as new 
stars are born, existing stars age and gas particles move 
in relation to the stars. This may drive additional non¬ 
equilibrium effects that we do not capture in our current 
simulations. 

Furthermore, since we have shown that UV radiation 
affects the global properties of galaxies, the inclusion of a 
fluctuating radiation field will also influence the galaxy di¬ 
rectly, in addition to any non-equilibrium effects that it may 
drive. Likewise, it will also be important to include fluctua¬ 
tions in metallicity due to chemical enrichment from stars, 
which are not included in the fixed-metallicity models that 
we present here. 
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